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PREFACE 


Our research on electromagnetic correlation techniques was initiated 
in 1965 to provide advanced experimental and analytical tools for investi- 
gation of turbulent or unstable flows. Tools to provide a level of scien- 
tific sophistication not possible in ad hoc engineering studies of launch 
vehicle configurations were sought. The application of electromagnetic 
correlation techniques has now been illustrated for the remote probing of 
sound source intensities in engine exhausts and for the detection of wind 
shears and turbulence. The above tasks require the coordination of the 
disciplines of sensor development, electromagnetic wave propagation, 
fluid mechanics, atmospheric science, physics, statistics, and data pro- 
cessing. Such a combination of disciplines was achieved through organiz- 
ing cooperative research programs with the Wave Propagation Laboratories 
of ESSA, the Department of Atmospheric Science of Colorado State University, 
the Department of Physics of the University of Oklahoma and MSFC's Compu- 
tation Laboratory. Cooperation in these programs meant a sharing of 
facilities and resources. 

This research review on electromagnetic correlation techniques covers 
contributions from inhouse engineers and scientists. Most of the^papers 
were contributed by undergraduate and graduate students during their 
training program at MSFC. This particular combination of cooperative 
research and cooperative training was initiated in 1966. 

It is hoped that this review will be stimulating to organizations 
which employ electromagnetic sensors for remote probing. Criticism of 
this review and discussions concerning individual papers with respective 
authors are invited and will be particularly helpful for the student 
trainees . 



W. K. Dahm 

Chief, Aerophysics Division 


H. J. Horn 
Deputy Director, 

Aero-As trodynamics Laboratory 
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CHAPTER I 


INFRARED MEASUREMENTS OF SOUND SOURCE INTENSITIES IN JETS 


by 

L. N. Wilson 

IIT Research Institute, Chicago, Illinois 

F. R. Krause and K. A. Kadrmas 
NASA - Marshall Space Flight Center, Alabama 


SUMMARY 


Direct measurements of local sound source intensities in turbulent 
jets are being attempted using an optical cross -correlation technique. 

It has been shown that optical cross-beam measurements can be related 
directly to the volume integral describing the source strength. This 
paper represents a study of this relationship. The application and 
limitations of the technique are discussed in detail, with reference 
to Lighthill's and, in particular, Ribner's mathematical formulations. 
The latter formulation is particularly suited to cross-beam measurements 
since it expresses source strengths in terms of volume integrals of mass 
density fluctuations. Area integrals of density fluctuations are being 
attempted using a new infrared cross-beam system. Progress to date on 
subsonic jets with this system shows promise of giving the desired dens- 
ity dependence. Measurements were made using the C0 2 absorption band at 
4.3j_i. This type of measurement should be applicable to estimations of 
Mach-wave emission. 


Symbol 


x 


NOMENCLATURE 

Definition 

position vector (x^; y^; z^) 
position vector (x; y; z) 

vector in direction of far-field observer (| x ; | 2 ; | 3 ) 



a 

p 

p 

m f 

k 

P 

P 

— > 

u 

a o 

a 

i 

I 

l 

T 

U C 


vector in direction of flow (£,r],£) 

mean density 

fluctuating density 

fluctuating mass fraction 

absorption 

mean pressure 

fluctuating pressure 

mean velocity vector (u; v; w) 

ambient speed of sound 

local speed of sound 

fluctuating signal 

mean signal 

path angle 

retarded time 

turbulence convection speed 


Subscripts 

A, B radiation beams of the crossed-beam system 

a contribution from absorption 

s contribution from scattering 


1-1. INTRODUCTION 


The mathematical formulation of the jet noise problem [1-1 , -2] has 
given us considerable insight into scaling laws associated with noise 
production, but little progress has been made in experimental verifica- 
tion of the predicted sound source strengths. The actual sound source 


2 



intensities as predicted by the theory are given in terms of turbulence 
flow properties which are extremely difficult to measure and even then 
must be incorporated into a complicated volume integral expression 
involving retarded time differences between statistical variables. 

Some attempts [1-3, -4], with a measure of success, have been made in 
estimating sound source strengths from hot wire measurements in low 
speed jet flows where the transformation to a moving frame of refer- 
ence allows retarded times to be ignored to a good degree of approxima- 
tion. At high speeds, this approximation is invalid; in fact, a 
singularity occurring at the sonic convection condition results in a 
prediction of infinite acoustic power radiated. A different type of 
noise source then dominates, and sound radiates from the jet in a manner 
analogous to Mach wave propagation [1-5]. 

An instrumentation technique, called the crossed-beam correlation 
method, has been recently developed [1-6] which automatically performs 
an area- integral of pertinent turbulent properties. This technique is 
similar in form to that required for a direct estimate of sound source 
intensity [1-7], and thus affords us the opportunity of measuring a 
part of the sound source volume integral directly. This paper covers 
the application of the technique to the Mach wave emission expected 
from supersonic jets and discusses progress on the experimental verifi- 
cation that has been achieved after the last review [1-8]. 


I- 2 . ANALYTICAL DEVELOPMENT 


A. Sound Source Intensities 

Ffowcs -Williams [1-5] has shown that, where turbulent jet convec- 
tion speeds are in excess of the ambient speed of sound 


( M c = — > !), 
o 

a pronounced finite sound peak occurs at cos 0 = 1/M C where 0 is the 
angle between the line to the observer and the flow direction (see fig- 
ure 1.1). Since this is the direction toward which Mach waves would ema- 
nate from supersonic bodies moving at a Mach number M c , the emission is 
analogous to the Mach wave emission from supersonically moving eddies, 
the latter behaving somewhat like solid bodies. The intensity of acous- 
tic radiation is then determined by the "shape" of the bodies (i.e., by 
the density gradient through the eddy). The sound source in this case 
is source-like (monopole); whereas, in subsonic flow the sound source is 
quadrupole in nature. The sound source strength is then determined 
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(through suitable mathematical manipulation) by density correlations in 
a fixed frame of reference, thus removing the complication resulting 
from measurements being required in the moving frame of reference as is 
the case of quadrupole noise. 

The far-field acoustic radiation from the shear layer of a turbu- 
lent jet is then predicted to be 


<P 2 > = 


1 

16n 2 


(X f - X) 2 (y f - y)2 

I — > — >1 6 

l X f “ X I 


{cWOy } 2 S(x) dx , 


where 


( 1 ) 




dp 

dt 


+ I, 


t 



( 2 ) 


S(x) is considered here to be the acoustic source strength per unit 
volume of turbulence. The other terms in the equations, apart from the 
mean velocity gradient, are dependent upon the relative position of the 
observer and the sound source. The velocity gradient can be measured 
independently of the source integral by other standard measurements. 
Ffowcs -Wil 1 iams writes the correlations in terms of the pressure gra- 
dients by assuming locally isotropic conditions: p = (l/a 2 )P. This 

technique requires a knowledge of the local speed of sound in the shear 
layer. If, however, we can measure density fluctuations directly, this 
approximation and requirement is not necessary. Ffowcs -Will iams also 
shows that the sound source integral is more properly written in a frame 
of reference aligned with and perpendicular to the Mach wave direction of 
propagation. 


B. Crossed-Beam Principle 

The operation of the cross ed-beam technique has been described 
elsewhere, and will therefore be discussed only briefly here. In 
principle, two narrow beams of selected wavelength radiation are posi- 
tioned such that they intersect at a point in the region of interest in 
the jet. If the beam intensity for each beam is modulated by a process 
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which is linearly related to the local gas density fluctuations, then 
the intensities can be expressed by 


i A = I A yVm f ) d h = 1 A f pdt] (3) 

i B = I B J (km f ) d^ = pd i, . (4) 


where i^, iB, IA and Ib are the fluctuating and mean signals, respec- 
tively, for each of the two beams A and B. If the two signals are 
cross -correlated over a sufficiently long integration time, the product 
mean value <iA ifi> will give information concerning only the fluctua- 
tions common to each beam (i.e., in the intersection region). In this 
way, local information on turbulence has been obtained without the 
necessity of inserting physical probes into the flow field [1-6]. From 
a more intuitive viewpoint (see figure 1.2), we might consider turbulence 
as consisting of discrete eddies that are slightly more or less trans- 
parent to the beams of radiation than the surroundings. As an eddy 
passes through a beam of radiation, the light intensity at the detector 
changes. When the two detector outputs are cross-correlated, only the 
information concerning eddies passing through the beam intersection 
point is retained. Signal fluctuations due to source instabilities, 
detector noise, etc., are normally uncorrelated, although these con- 
tributions do affect the overall signal-to-noise ratio of the system 

[1-9]. 


Cross -correlation of equations (3) and (4) gives 


G(x) = 


<1 A 1 B > 

<I A ><I B > 


fj <p(x,y+T],z,t) p(x+|,y,z,t+T) >d T1 dg, (5) 


where i A and ig are correlated with an arbitrary time delay x between 
signals and | is a displacement of beam B relative to beam A in the 
direction of U c . This equation expresses the correlation between two 
line integrals. If we can write equation (5) as 


G(x) = 


<p(x,y,z,t) p(x+f* ,y+T), z+£ , t+f)> dq d(; 


( 6 ) 
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or 


G(x) = <p(x,y,z,t) J J <p(x+|,y+T), z+^, t+x) dq df;> 
point area 


(7) 


then, we are dealing with the correlation between the area integral of 
points in a plane (x + | = constant) with an outside point (x,y,z). 

Such a cross-correlation has been assigned the term "point-area correla- 
tion," and is the type occurring in equation (1). 

If the turbulent flow is everywhere homogeneous, then the point- 
area formulation is obviously proper. If departures from homogeneity 
are pure odd functions of space coordinates about the point (x,y,z), 
then these departures tend to integrate out and the point-area concept 
is still valid. In any event, the final proof of the point-area correla- 
tion assumption can come only from experiment. Two tests for this proof, 
discussed in detail in reference 1-8, are summarized here. 

(1) The spatial resolution of the measurements across stream- 
lines should be limited only by the outside "point" which, in practice, 
is the optical beam diameter. 

(2) Since the cross-correlation of signals i^ and ig will 
correspond only to contributions from turbulent fluctuations normal to 
the beam plane, the energy spectrum should approximate the three-dimen- 
sional spectrum of turbulence. This is in contrast to a two-point hot 
wire measurement of the so-called one-dimensional spectrum which contains 
contributions due to cross-flow components. As a result, the energy 
spectrum for crossed-beam data should approach a -11/3 variation in the 
high frequency portion. 

Therefore, assuming the validity of the point-area correlation 
approximation, equation (7) in abbreviated fashion becomes 


G(x) 



p(x + £, t + t)> d rj d£. 


( 8 ) 
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C. Application to Sound Source Intensity Measurements 


The similarity between equation (1) and equation (8) is striking. 

It would appear that only one further integration (point-area to point- 
volume) is required to obtain the sound source intensity. There are, 
however, two main points to be considered carefully in making the com- 
parison. They are as follows: 

(1) The radiation extinction process must be directly 
gas-density dependent. 

(2) Time derivative and not space derivative correla- 
tions are required, and the correct retarded time 
difference must be used. 

We will now consider these two points in greater detail. 

1 . Extinction Process 

Several extinction processes have been used in crossed-beam 
studies in the past. These include scattering from natural and arti- 
ficial tracers and the use of the absorption of air in the ultra-violet 
portion of the spectrum (Schumann-Runge bands of oxygen). Each of these 
has proven to be unacceptable for sound source measurements although 
they are useful in other turbulence parameter studies. More recently, 
instrumentation has been developed to measure absorption by naturally 
occurring C0 2 in air using the near infrared portion of spectrum. The 
fundamental vibration band of C0 2 near 4.3 microns is the most promising 
due to the absence of water vapor absorption in this region. Only 
naturally occurring C0 2 can be used, since only then will the optical 
fluctuations measured be truly representative of the gas density fluc- 
tuations in the jet. The use of the infrared spectrum also reduces the 
contributions to the signal fluctuations from ever present particle 
contaminants in the flow. 

The extinction process in general can be written as 


i + I = I D exp -/ (p + p )di 


= I exp 
o 


di 


exp 


p d£ 


= I exp 


p d^ 


I 
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Now, for small perturbations 


exp 



p 




1 



p du. 


Thus , 


i = 


-I 



p 


This was the form presented in a previous section. 

It should be repeated at this point that p is analogous to 
km-g. This means that a comprehensive study must be made using a unique 
absorption calibration cell to determine the dependence of k on the 
static and dynamic pressure and static temperature. Therefore, knowing 
Sk/SP and Sk/ST, the experimenter can choose the optimum wavelength and 
bandpass settings in the absorption band to discriminate against the 
temperature dependence of k. These experiments can be made only in a 
controlled environment free of a contribution due to scattering to the 
observed signal. 

In principle then, it is feasible to use infrared absorption 
by naturally occurring C0 2 for the measurement of gas-density fluctua- 
tions. The preliminary results are given in the next section. 


2 . Time Derivative Requirements and Retarded Time Effects 

The point-area correlation in equation (6) must be in terms of 
the time derivatives of the gas-density fluctuations in order to obtain 
an experimentally accessible form of equations (1). 

For stationary time series, the correlation of time derivatives 
is identical with the time lag differentiation of the correlation. 


SpC?, t) Sp(x+£, t+r) 

" < — — > = 


St 


St 


S 2 <p(x,t) p(x+|, t+i)> 

__ 


( 9 ) 


s 2 g (x) _ r r s 2 < P (x,t) p(x+^, t+i)> 

St 2 J J St 2 
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( 11 ) 


S(x) 


ci 2 G (?) 

St 2 


dll- 


We can also simplify the retarded times in equation (1) by the approxi 
mation 


— > 
x 


?! 



( 12 ) 


The order of time lag differentiation and space integration (equation 
(11)) may be interchanged in the special coordinate system g = ({■]_; | 2 ; 
£ 3 ) where the planes = constant are parallel to the Mach wave fronts 
that are generated at the source point x*, the reason being that retarded 
times are constant along such planes within the far-field approximation. 
This type of crossed-beam orientation is possible, but does not allevi- 
ate the problem since 8 is not known a priori . 

Now, the time delay for the orientation in figure 1.2 is 


T 






This time delay is also precisely that which allows an eddy to move 
from the upstream beam to the downstream beam along a streamline in a 
reference frame oriented in the flow direction | = ((■; q; £). In other 
words, the retarded time is correct for the area near the "common" 
streamline, but away from the streamline r is a function of the particu- 
lar points t| and In the first approximation, the contributions should 

cancel as in the case for cross-flow components. Therefore, we can write 
equations (10) and (11) as 


<3 2 G (?) P P a 2 <p(?,t) o(x+£, t+x) 

St 2 J j 0 t 2 


(13) 


S(x) 



c> 2 G (?) 
St 2 




(14) 
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and 



— > 


— > 

— > 

x f 

X 


x f 

- x - 1 


( 15 ) 


Now, S(x), being an experimentally accessible quantity under the above 
assumptions, sound source intensities can be measured. 


1-3. EXPERIMENTAL STUDIES 


The first tests* conducted with the infrared crossed-beam system 
were made in the shear layer at an axial location of X/D = 4 of a sub- 
sonic jet so that comparisons could be made with more standard measuring 
techniques. A photograph of the one-inch diameter subsonic jet facility 
and infrared crossed-beam system is shown in figure 1.3. 

The radiation sources currently in use are glowbars operating at a 
temperature of 1100°K (equivalent blackbody temperature). Cass egrainian 
optics are used to focus the radiation from the sources into the jet 
flow field with an image dimension of 2mm in the flow direction and 
1 mm in the transverse direction. The radiation is then collected, using 
similar optics and refocused on the entrance slit of a McPherson Model 
218 monochromator. The experimenter can then choose the center wavelength 
and bandpass settings desired in the 4.3 micron absorption band of C0 2 . 

The bandwidth initially used was typically 0.08 microns. The radiation 
level is monitored at the exit slit of the monochromator by using LN 2 
cooled Indium-Antimonide detectors. The fluctuating signals obtained 
are amplified and recorded for a detailed analysis. Simultaneously, the 
signals are reduced on-line using a PAR correlator. A computer program 
is available to reduce the digitized data tapes to obtain the correla- 
tions and spectra. 

The cross-correlations (zero beam separation) and mean levels 
obtained while traversing the spectral region from 4.1 to 4.5 microns 
are shown in figure 1.4. Both the mean level traverse and cross-corre- 
lation traverse follow the shape of the C0 2 absorption curve as expected, 
but cross -correlations were recorded outside the absorption band where 
the absorption should fall to zero. It is felt that this is due to the 
contribution from scattering to the total signal since the jet has con- 
tamination from oil mist due to the compressors. In this case, the 


The data presented here were kindly supplied by Dr. R. J. Damkevala of 
IIT Research Institute. 
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signal is obviously not entirely due to C0 2 absorption, and therefore 
some method must be devised to treat this type of problem. 

If we consider the fluctuating signals of beam A and beam B to be 
composed of contributions due to absorption and scattering, then 


i A (A, . ..) = i A (A, . ..) + i A (A l5 ...) 


and 


i-rj (A, . . . ) ~ ir> (A, • . • ) + i-D (A x , . . . ) , 
o u a D s 


(16) 


where we have assumed that the scattering contribution is a linear con- 
tinuum across the absorption band and can thus be approximated at A^ 

The cross-correlation of the two composite signals yields 

^ i A ( A, ...) ig ( A, • • • ) A* — r (A, » • . ) C A 9 ...) 

+ < (A, . . . ) ig (A l5 . . . ) > + < i^ (Ai, . • • ) ig (A, . . . )> 

Si S S 3 

+ < i^ g (A-|_, . . . ) ig s (A ls ...)>. (17) 

Therefore, to obtain 

^ 1a (A, * • • ) tg (Aj • • • ) A’ , 
a D a 

we must know the three remaining terms on the RHS of equation (17). 

By definition, 


Qaa< A ’ A > 


<i-Ao • • * ) i-Ba • . • ) 

3 3 

< X A > < I B > 


(18) 
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and 


Q ss (Al ’ Al) 


<iA s (Ai, ...) 1 b s (Ai, ...)> 
< I A > < S > 


(19) 


Also, 


Q as = (A, A,) 


<iA a (A, . ..) iB s (A ls 


..)> 


< X A > < X B > 


< i A 0 (Ai, . ••) iB_(A, ...)> 



( 20 ) 


since both processes are governed by the turbulence convection and are 
thus well correlated. Now, setting beam A at }\- L and beam B at A, cross- 
correlating the two signals, we get 


^ (Aj_, •••) ig (A, ...)A > 

< X A > < T B > 


= Q as ^ A ’ A i^ + Q S s ( A l’ A i) 


(21) 


or 


Qas< A ’ A ^ “ 


<\(Ai, ...) ig(A, ...)> 

< I A > < > 


- Q ss (Ai»A 1 ). 


( 22 ) 


Finally, we can write 


q aa a.» = 


i* (A, • • • ) i-n (A, •••)!> 

<A><V> 


- 2 


<i^ ( A^j • • • ) ig (A, . . . )> 
< i A > C T B > 


+ Qgg (Al, Aj) , 

which is the contribution due to absorption only (see figure 1.4). 


(23) 
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Figure 1.5 shows a radial profile measured at a center wavelength 
of A = 4.24 microns. The resolution of the crossed-beam system may be 
assessed by comparing the width of this intensity peak with the shear 
layer width. The consistent variation of the two-beam product mean 
value is evidence of resolution capabilities, and supports the earlier 
contention of the validity of the point-area correlation concept. 

Fourier transforms performed on the cross-correlations are pre- 
sented in figure 1.6. Once again the point-area concept appears to be 
supported by experimental evidence, since the high frequency portion of 
the spectrum does approximate the -11/3 power law very well. 


1-4. CONCLUSIONS AND FUTURE PLANS 


The very promising first results presented in this paper have given 
confidence that the crossed-beam technique will allow measurements of 
sound source intensities to be made in jet flows. Toward that end, 
experiments are being initiated to measure two-beam correlations in a 
heated subsonic jet, and measurements have already begun in a specially 
designed absorption calibration cell to verify analytical predictions, 
with regard to understanding pressure and temperature effects. These 
measurements are to be made by varying the absorber mass, static and 
dynamic pressure and static temperature independently. 

With this information and the already proven ability to obtain 
remotely point-area correlations that are not accessible to probe 
instruments, the crossed-beam system will have the potential of locating 
and measuring, for the first time, the turbulent sound source intensities 
in jet flows. Such measurements, taken without disturbing the flow, will 
not only lend support to analytical predictions, but also aid in the 
design and evaluation of noise suppression devices. 
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FIGURE 1.1. EXTENSION OF ACOUSTIC WAVES INTO FLOW 




FIGURE 1 . 2 . CROSS-BEAM APPROXIMATIONS OF LOCAL SOUND SOURCE INTENSITIES 





FIGURE 1.4. DECOUPLING OF ABSORPTION AND SCATTERING 





2-16 m/sec 



FIGURE 1.5. RADIAL PROFILE OF CROSSED BEAM INTENSITY 
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CHAPTER II 


OPTICAL DETECTION OF MOVING STRIATIONS 
WITH THE CROSS -CORRELATION TECHNIQUE 


by 

J. B. Stephens and J. B. Thomison, Jr. 
NASA - Marshall Space Flight Center, Alabama 


SUMMARY 


Optical observations of emission variations in a plasma generated 
by a glow discharge have been achieved using the cross-beam correlation 
technique. Classification of these emission variations at a point into 
periodic, narrow band, and random processes, is discussed, along with 
the identification of convection speeds and feedback loops of the system. 


II- 1. INTRODUCTION 


The experiments documented in this paper were designed to develop 
cross-beam methods for studying local photo-chemical processes in a 
plasma. Previous investigations of glow discharges by Abria [II-l] and 
later by others [II-2] indicated that it is composed of moving optical 
fluctuations. More recent photographic experiments by Garscadden and 
Lee [II-3] show that these fluctuations can move both toward the cathode 
and toward the anode. Our tests differ from previous investigations by 
applying the cross-beam correlation technique which was developed in 
the wind tunnel by Fisher and Krause [II-4] and in the atmosphere by 
Stephens et al. [ II - 5 ] to measure wind speeds. Velocity and frequency 
information has been obtained by extending the cross-beam technique from 
extinction processes to emission processes. 
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II- 2. PHOTOGRAPHIC EVALUATION OF GLOW DISCHARGE PHENOMENA 


The plasma instabilities that produced the optical fluctuations are 
defined to be striations in the positive column of the glow discharge 
(see Fig. II. 1). Stationary striations, which are packets of ionized 
gas separated by dark spots, normally occur at tube pressures below 
1 mm Hg (see Fig. II. lb). Figure II. la is an example of moving stria- 
tions. Although this picture was taken at 1/100 of a second, fluctua- 
tions are not obvious, and the glow appears as a continuum. Figure II. 2 
shows the test setup as it appeared immediately before a high-speed 
movie was taken. 

If a striation contains an excess of positively ionized gas, it 
should move toward the cathode. Conversely, if a striation contains an 
excess of electrons, it should move toward the anode. This is in accord 
with the current views of plasma physics [ I I- 3 ] . These two cases are 
illustrated in Figures II. 3 and II. 4 taken at 1/400 of a second. In 
Figure II. 3 the smaller striation moves primarily toward the cathode, 
reverses its direction for a short time, and moves backward toward the 
anode. In Figure II. 4 the movement of the striation tends entirely 
toward the anode. This is in keeping with the results of Garscadden 
and Lee [II-3] which show existence for forward- and backward -moving 
striations . 

High-speed movies were taken of this phenomenon, but, unfortunately, 
movies provide only qualitative descriptions since the photographic film 
is a logarithmic detector. Therefore, one cannot see a fine structure in 
the photographs of the large semi-stationary ionization. The present 
state-of-the-art in photography restricts measurements to those occurring 
below 1000 Hz. For example, little or no detailed information can be 
obtained about a glow discharge in argon since its frequency range is 
from 1000 Hz to at least 10,000 Hz. 


II- 3. CROSS-BEAM EXPERIMENTS 


Both qualitative and quantitative results have been achieved by 
using two crossed photomultipliers. These photomultipliers are linear 
over their entire range and have a frequency response better than our 
electronics . 

The experimental setup is shown in Figure II. 5. Our measurements, 
which were restricted to the positive colunuij were made at an optical 
wave length of 3100 A with a bandpass of 50 A. The focal length of the 
optics was 11.3 cm with a diameter of the beam at intersection of 1 mm. 
The work was done in a discharge tube of a diameter of 9 cm and with 



stainless steel electrodes 82 cm apart. The major portion of our effort 
was concentrated in Argon at the pressure range from 1 mm Hg to 25 mm Hg, 
with some work also being done in C0 2 and air. 

The outputs of the optical collectors were then fed into the 
Princeton Applied Research (PAR) analog correlator to obtain a covar- 
iance reading. This on-line analog computer assumes that no D.C. shifts, 
trends or gain instabilities were introduced by the electronics. 

Further, the covariance thus obtained will become uncertain as soon as 
the correlation drops below 10 percent of the mean signal power. 

Within the above assumptions, the analog computations of the PAR 
correlator should form an estimate of the temporal covariance function, 
which is defined as follows: 


T 

“ (t) = -~J (X(t) - X)(Y(t- x ) - Y) dt, 
o 

where 

x(t) = instantaneous fluctuations of intensity from the 
undelayed channel 

Y(t - t) = instantaneous fluctuations from the delayed channel 
X and Y= average values of the fluctuations 
t = delay time set during data reduction 

t = real time 

T = total integration time 

N = normalization factor. 

This covariance enables one to obtain the best estimate for the 
amount of time required for the striation to traverse the transit 
distance of the two beams (see Fig. II. 5). The covariance also provides 
a systematic approach to the classification of narrow- and broad-band 
processes while systematically treating an entire ensemble. 

As shown in Figure II. 6 several types of wave form were encountered. 
The first signal is easy to classify since it is obviously periodic, but 
the remaining signals require some kind of statistical description. 

Figure II . 7 shows how the top three signals can be classified by use of 
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an autocorrelation curve. When the maximum amplitude of the peaks does 
not decay as a function of time delay, then the process is periodic. 

When the peak decay is a function of time delay the process is either 
narrow-band or broad-band, depending on the rate of correlation decay. 
The broad-band process occurs where there is only a fundamental peak 
and the remaining peaks are washed out. In the narrow band process 
the higher order peaks decrease in amplitude as the order increases. 
Physically, these three processes tell you that; 

(1) If the signal is periodic, then the signal time 
history is recurring at exact intervals of time. 

(2) If the signal is narrow-band, then the time history 
will repeat itself for a finite number of periods 
until it is damped out. 

(3) If the signal is broad-band, then the signal time 
history will not repeat itself. 

Thus, the kind of process provides a means of classification of the 
internal consistency of the structure. 

The fourth signal on Figure II. 6 affords a more challenging test 
for the correlation technique as shown in Figure II. 8. It was found 
that by filtering the signal it could be subdivided into a high fre- 
quency and a low frequency component. The low frequency component 
was due to periodic fluctuations and moved in the direction of the 
cathode. The high frequency component was a narrow-band process and 
had a velocity in the direction of the anode. This could well be 
thought of as an ion striation (moving toward the cathode) with a sub- 
structure composed of electron waves (moving toward the anode). The 
direction of flow was determined by starting with intersecting beams 
and increasing the beam separation. As the separation increased, the 
transit time increased positively or negatively according to the flow 
direction. 

Another kind of measurement performed was the identification of the 
feedback loop between the fluctuations of the discharge current and the 
optical fluctuations in emission. This was done by cross -correlating 
the optical fluctuations with the current fluctuations. When the beam 
separation increased, the transit time increased (bottom correlation 
curve in Figure II. 8). Thus, it was possible to close the feedback loop 
and say that the fluctuations in the current were due to the flow of 
electrons in the discharge. 

When one changes the potential across the discharge, the discharge 
is unstable until it reaches thermal equilibrium. All of our measure- 
ments were made after a short period of time had elapsed, and it was 
confirmed that all parameters were stable. 
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II-4. PERIODIC STRIATIONS 


The frequencies of the periodic fluctuations ranged from 10 Hz to 
50 KHz and the velocities from 0 to 30 km/sec. The results are based 
on over 200 measurements, which are summarized in Figures II. 9, II. 10, 
II. 11 and Table II. 1. 

During the course of our measurements, it was noticed that a 
periodic waveform would appear out of what had been a random narrow- 
band process. This periodic waveform was stable for the conditions 
under which it was observed, but any slight change in conditions 
resulted in a reversion to the random narrow-band process. The 
periodic process was observed in Argon over a wide range (1 mm - 20 mm) 
of pressure, but only certain discharge voltages for each pressure. It 
was also noticed that the amplitude of the periodic waveform was much 
larger than that for the associated random narrow-band process. 

The wave length of the periodic waveform seemed to fall approxi- 
mately into two categories (see Fig. II. 10): 

(1) Wave length equal to length of tube; 

(2) Wave length equal to one half of tube length. 

The velocity and frequency of the periodic striations are prac- 
tically independent from the tube pressure within the range between 1 
and 20 mm Hg as illustrated in Figures II. 9 and II. 10. The kind of 
process, the kind of striation (ion or electron), and the gas composi- 
tion all have a first order effect on the velocity, frequency, and 
wave length of moving striations. Even though we could not clearly 
establish a pattern, there appeared to be times during which only one 
type of striation was present. Ion and electron striations did not 
always occur simultaneously. 


II-5 . STANDING STRIATIONS 


The second distinct class of striation that was observed in the 
glow discharge is a standing striation. We always found that the fluc- 
tuations in emission associated with standing striations were broad band. 
The first zero crossing of the auto-correlation function was taken as an 
estimate of the time period associated with the center frequency of the 
broad band process. This period was divided into the measured velocity 
to obtain a characteristic wave length. 
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Although the striations appear standing to the observer, the cross 
correlation tests indicated that they have small velocities about a 
factor of a hundred less than those associated with the periodic stria- 
tions. The direction of these small velocities seemed to follow only 
the electron flux. 


II-6. CONCLUSIONS AND RECOMMENDATIONS 


Since our results for the velocities and frequencies of moving 
striations are in keeping with those to be found in the literature, we 
have confidence in the cross-beam technique to retrieve this information 
in emission. 

This research has demonstrated the feasibility of the cross-beam 
technique for successful retrieval in emission of (1) phase and group 
velocities along the main current, (2) separation of periodic and random 
narrow-band components, (3) classification of optical disturbances using 
length and time scales in moving frame, (4) classification of flow 
instabilities through statistical evaluation of changes in surroundings, 
and (5) correlation between feedback loops and the control system (e.g. , 
correlation of fluctuations in emission and fluctuations in current). 
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FIGURE II. 2. EQUIPMENT SETUP 
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FIGURE II. 3. ION STRIATION 
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FIGURE II. 7. CORRELATION OF SINGLE PHOTOMETER OUTPUTS FOR TYPICAL FLUCTUATIONS 




ELECTRODE [ 


ELECTRON FLOW 


0 I ELECTRODE 
CURRENT 


< 8 > ® 


ION FLOW 



FIGURE II. 8. SEPARATION OF ELECTRON AND ION STRIATIONS 
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FIGURE II. 9. EFFECT OF STRIATION DIRECTION ON ITS VELOCITY 
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WAVE LENGTH (METERS) 





FIGURE II. 10. QUANTUM EFFECT OF WAVELENGTH 
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TABLE II. 1. EXPERIMENTAL DATA 


GLOW DISCHARGE IN ARGON 


RUN NUMBER 

» 20.0 

PRESSURE 

. 5 

;o 





POWER 

RESIST 

PREQ 

VEL 

WLG 

DlS I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHz 

KM/SEC 

meters 

RA 

VOLTS 

USEC 

USEC 

8.100 

36*000 

17*857 

12*097 

*677 

15.0 

540*C 

56.0 

24.8 

33.605 

7.954 

18*519 

12*295 

*664 

65.0 

517*0 

54.0 

24*4 

59.800 

4.522 

19*231 

12.500 

*650 

U5. 0 

520 • C 

52*0 

24.0 

90.255 

3.315 

20*833 

12*500 

•600 

185.0 

547*0 

48.0 

24.0 

125.560 

2.716 

20*833 

13*636 

* 655 

215.0 

584 . c 

48.0 

22.0 

166.950 

2*377 

22*222 

13*889 

*625 

265.0 

630»C 

45.0 

21.6 

211.050 

2*127 

23*810 

15.000 

• 630 

315.0 

670*0 

42.0 

20.0 

264.625 

1*986 

26*316 

15.789 

•600 

365.0 

725*0 

38.0 

1.9.0 

321.625 

1*867 

27*027 

21*429 

• 793 

415.0 

775*0 

37*0 

14.0 


RUN NUMBER 

* 23*0 

PRESSURE 

* 4 

.8 





POWER 

RESIST 

FREq 

vel 

WLG 

DIS I 

DIS V 

PER 

TAU 

WATTS K 

OHMS 

KHz 

KM/SEC 

meters 

RA 

VOLTS 

USEC 

USEC 

25 . 650 

10*260 

1*600 

*498 

*311 

50.0 

513*0 

625.0 

603.0 

51.500 

5*150 

2*353 

•744 

*316 

100.0 

515*C 

425.0 

403.0 

82.500 

3*667 

2*857 

•990 

*347 

150*0 

55o*0 

350.0 

303.0 

117.000 

2*925 

25*641 

15*000 

.585 

200*0 

585 • 0 

39.0 

20.0 

155.000 

2*480 

25*641 

15.385 

*600 

250.0 

620*0 

39.0 

19.5 

297.600 

1*860 

28*571 

18*182 

*636 

400.0 

744*0 

35.0 

16.5 

351 . 450 

1*736 

25*641 

18.182 

• 709 

450.0 

781*0 

39.0 

16.5 

298.000 

1*862 

45 * ^55 

15.957 

*351 

400.0 

745*0 

22.0 

18.8 

247.450 

2*020 

23*256 

16*667 

•717 

350.0 

707*0 

43.0 

18.0 

200.400 

2.227 

22*472 

16.216 

*722 

300.0 

668*0 

44.5 

18.5 

158.000 

2*528 

21*598 

15*957 

• 739 

250*0 

632.Q 

46.3 

18.8 

118,200 

2*955 

20*964 

15.385 

*734 

200*0 

59l *o 

47.7 

19.5 


RUN NUMBER 

« 24.0 

PRESSURE 

■ 4 

.9 





POWER 

RESIST 

FREQ 

VEL 

WLG 

DIS I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHz 

KM/seC 

meters. 

HA 

VOLTS 

USEC 

USEC 

30.900 

12*360 

22*124 

16.760 

.758 

50.0 

618*0 

45.2 

17.9 

61.800 

6*180 

25*000 

17.143 

.686 

100.0 

618*0 

40.0 

17.5 

95.550 

4*247 

25*510 

18.072 

.708 

150.0 

637*0 

39.2 

16.6 

134.000 

3.350 

26 . 1 78 

18.750 

.716 

200.0 

67o*0 

38.2 

16.0 

175.750 

2*812 

26*667 

18.750 

•703 

250.0 

703*C 

37.5 

16.0 

221.700 

2*463 

28*090 

17.544 

• 625 

3$0 • 0 

739*0 

35 . 6 

17*1 

269.850 

2*203 

27*174 

18.750 

• 690 

360.0 

771*0 

36.8 

16.0 

321.200 

2*007 

27*397 

18.987 

• 693 

400.0 

803*0 

36.5 

15.8 

376.650 

1*860 

27*174 

20*408 

•751 

460.0 

837*0 

36.8 

14.7 

440.000 

1*760 

29.240 

21.127 

• 723 

500.0 

880*0 

34.2 

14.2 
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TABLE II. 1. (Continued) 


RUN NUMBER 

= 25.0 

PRESSURE 

: * 9. 

.8 





POWER 

RESIST 

FREQ 

VEL 

WLG 

DIS 1 

DIS V 

PER 

TAU 

WATTS 1 

< OHMS 

KHz 

KM/SEC 

meters 

MA 

VOLTS 

USEC 

USEC 

8.970 

39.867 

14*184 

10*526 

•742 

15.0 

598 • C 

70*5 

28.5 

29.500 

11.800 

25*510 

11*321 

*444 

60.0 

590 • C 

39.2 

26.5 

59.300 

5.930 

26*385 

11*673 

*442 

100.0 

593*0 

37.9 

25.7 

90.750 

4.033 

22*779 

12*146 

• 533 

150.0 

605*0 

43.9 

24.7 

125.600 

3-140 

24.510 

13*453 

*549 

200.0 

628*0 

40*8 

22.3 

160.000 

2*560 

35.088 

13*514 

•385 

250.0 

640*0 

28.5 

22.2 

205.500 

2.283 

56.497 

13.393 

•237 

300.0 

685*C 

17.7 

22.4 

253.050 

2*066 

59-524 

14.778 

•248 

350.0 

723*0 

16.8 

20.3 

303.600 

1*897 

60*241 

15.789 

•262 

400.0 

759*0 

16.6 

19.0 

358.200 

1*769 

65*789 

16-575 

• 252 

450.0 

796*c 

15.2 

18.1 

417.500 

1*670 

76-336 

16*760 

•220 

500.0 

835*0 

13.1 

17.9 

139.200 

3.480 

54.348 

15.707 

•289 

200.0 

696*0 

18.4 

19.1 

101.400 

4.507 

60* 606 

15.228 

• 251 

150.0 

676*0 

16.5 

19.7 

66.500 

6.650 

57*143 

14.493 

•254 

100.0 

665*0 

17.5 

20.7 

32.700 

13.080 

21*739 

15.464 

• 711 

50.0 

654 • C 

46.0 

19.4 

RUN NUMBER 

power 

* 28.0 
RESIST 

PRESSURE • 5 

FREQ VEL 

a 

WLG 

D IS I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHz 

KM/SEC 

meters 

MA 

VOLTS 

USEC 

USEC 

30-400 

12*160 

24*752 

16.854 

• 681 

50.0 

608*0 

40*4 

17.8 

60.800 

6*080 

25-000 

16.484 

• 659 

100.0 

608*0 

40.0 

18.2 

94.800 

4.213 

25.316 

16.216 

• 641 

160.0 

632*c 

39.5 

18.5 

131.400 

3-285 

25.000 

16.375 

• 663 

200.0 

657*0 

40*0 

18.1 

171.750 

2.748 

26*667 

17.442 

.654 

250*0 

687*0 

37.5 

17.2 

215.400 

2.393 

26.525 

18.073 

.681 

300.0 

718*0 

37.7 

16*6 

262.150 

2.140 

26 . 596 

18.519 

.696 

350.0 

749*0 

37.6 

16.5 

274.000 

1.712 

27.174 

18.519 

• 681 

400.0 

685*0 

36.8 

16.2 

367 . 650 

1*816 

28.409 

18.987 

• 668 

460*0 

817*0 

35.2 

15.8 

428.000 

1.712 

28.736 

20.134 

.701 

500.0 

856*0 

34.8 

14.9 

RUN NUMBER 

power 

« 29.0 

resist 

PRESSURE * 2 

FREQ VEL 

.0 

WLG 

DIS I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHz 

KM/SEC 

meters 

RA 

VOLTS 

USEC 

USEC 

23.250 

9*300 

20 * 408 

15.306 

• 750 

50.0 

465*0 

49.0 

19.6 

53.800 

5*380 

21.008 

14.925 

• 710 

100.0 

538*0 

47.6 

20*1 

91 >050 

4*047 

21*930 

14.778 

*674 

150.0 

607*0 

45.6 

20.3 

134.600 

3.365 

25.126 

14.851 

• 591 

200.0 

673*0 

39.8 

20.2 

184.250 

2*948 

24*213 

14.354 

• 593 

250.0 

737*0 

41.3 

20.9 

239.400 

2.660 

22*989 

15.789 

• 687 

300.0 

798*0 

43.5 

19.0 

303.800 

2.480 

25*316 

16.854 

• 666 

350.0 

868*0 

39.5 

17.8 

370.400 

2.315 

26 • 882 

17.964 

• 668 

400.0 

926*0 

37.2 

16.7 

449.550 

2.220 

28*571 

19*737 

• 691 

460.0 

999*0 

35.0 

15,2 

531.500 

2.126 

29.412 

20.000 

• 680 

500.0 

1063*0 

34.0 

15.0 

RUN NUMBER 

power 

» 30.0 

resist 

PRESSURE <■ 9 

FREQ VEL 

.8 

WLG 

DIB I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHZ 

KM/SEC 

meters 

RA 

VOLTS 

UBEC 

USEC 

129.584 

3.912 

23*041 

15.957 

• 693 

182.0 

712*0 

43.4 

18.8 

141.400 

3.535 

22-831 

16.854 

.738 

200.0 

707*0 

43.8 

17.8 

105.900 

4.707 

22*989 

16.575 

• 721 

150.0 

706*0 

43.5 

18.1 

70.900 

7*090 

21*598 

16.129 

.747 

100.0 

709*0 

46.3 

18.6 

54. 450 

9.680 

22*371 

15.625 

• 698 

75.0 

726*0 

44.7 

19.2 

101.700 

4.520 

22*523 

16.216 

.720 

150.0 

678*0 

44.4 

18.5 

139.600 

3.490 

25*253 

16 . 667 

. 660 

200.0 

698*0 

39.6 

18.0 

178.Q00 

2.848 

24*938 

17-341 

• 695 

250.0 

712*0 

40.1 

17.3 

222*600 

2.473 

25-840 

17.857 

• 691 

300.0 

742*0 

38.7 

16.8 

271.250 

2.214 

26*667 

18.750 

*703 

350.0 

775*0 

37.5 

16.0 

325.200 

2.032 

28*169 

20 .Q 00 

-•710 

400.0 

813*0 

35.5 

15.0 

380.250 

1*878 

29*851 

20.408 

*684 

450.0 

845*0 

33.5 

14.7 

440.500 

1*762 

30*120 

21.127 

• 701 

500.0 

881*0 

33.2 

14,2 

385 . 650 

1*904 

30*303 

21*127 

.697 

450.0 

857*0 

33.0 

14.2 
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TABLE II. 1. (CONTINUED) 


RUN NUMBER 

a 26*0 

PRESSURE 

2 

*0 





power 

RESIST 

FREQ 

VEL 

WLG 

DlS I 

DIS V 

PER 

TAU 

WATTS 

K OHMS 

KHz 

KM/SEC 

meters 

MA 

VOLTS 

USEC 

USEC 

41.905 

5«8oO 

13»605 

10.169 

.747 

85.0 

493*0 

73.5 

29.5 

51.100 

5.110 

14.599 

10,791 

.739 

100. 0 

511*0 

68.5 

27.8 

86*250 

3.833 

15*385 

11*628 

• 756 

150.0 

575*0 

65.0 

25.8 

128 . 600 

3.215 

17*094 

12.500 

*731 

2Q0.0 

643*0 

58.5 

24.0 

176.250 

2*820 

18*349 

13*453 

.733 

250.0 

705.0 

54.5 

22.3 

231.000 

2.567 

21*008 

14*634 

*697 

300.0 

770*0 

47.6 

20,6 

297.500 

2.429 

23*697 

16.575 

.699 

35Q.0 

850*0 

42,2 

18.1 

374.000 

2.337 

26,667 

18.519 

• 694 

400.0 

935*0 

37.5 

16.2 

463.500 

2.289 

6*849 

2.222 

*324 

450.0 

1030*0 

146.0 

135.0 

553.000 

2.212 

7.634 

2.655 

• 348 

500.0 

1106*0 

131*0 

113.0 

397.200 

2*482 

6*557 

2.344 

• 357 

400.0 

993*0 

152.5 

128.0 

337 . o5o 

2.751 

5*650 

'2. 027 

-.359 

360.0 

963*0 

177.0 

-148.0 

146.Q00 

3.650 

27.634 

16.043 

.581 

200.0 

730*0 

36.2 

18.7 

195.500 

3.128 

26.247 

15.625 

• 595 

250.0 

782*0 

38.1 

19.2 

252.000 

2*800 

25*000 

16.484 

• 659 

300.0 

840*0 

40.0 

18.2 

311.500 

2.543 

15.674 

21.739 

1.387 

350.0 

890*0 

63.8 

13*8 


RUN NUMBER 

a 27.0 

pressure 

B 

1 

.9 







power 

resist 

FREQ 

VEL 

WLG 

DIS I 


DIS V 

PER 


TAU 

WATTS 

K OHMS 

KHZ 

KM/SEC 

METERS 

MA 


VOLTS 

USEC 

USEC 

163 

59o 

3*710 

33 

333 

20 

.690 

• 621 

210 

0 

779*0 

30 

.0 

14.5 

99 

000 

4.400 

32 

258 

20 

•000 

• 620 

150 

0 

660*0 

31 

.0 

15,0 

62 

100 

6.210 

29 

851 

20 

•000 

• 670 

100 

0 

621. C 

33 

.5 

15,0 

26 

200 

10*480 

27 

027 

15 

•625 

,578 

50 

0 

524*0 

37 

*0 

19.2 

1 

207 

1*932 

25 

000 

15 

• 625 

• 625 

35 

0 

48*3 

40 

.0 

19.2 

25 

000 

10*000 

24 

938 

16 

• 667 

.668 

60 

0 

500*0 

40 

*1 

18*0 

55 

300 

5*530 

25 

510 

16 

• 667 

• 653 

100 

0 

553*0 

39 

• 2 

18.0 

93 

150 

4*140 

26 

810 

16 

•216 

• 605 

1B0 

0 

621*0 

37 

• 3 

18.5 

137 

600 

3*440 

27 

548 

16 

• 484 

• 598 

200 

0 

688*0 

36 

*3 

18,2 

189 

000 

3*024 

29 

240 

16 

*667 

•570 

250 

0 

756*0 

34 

*2 

18.0 

249 

600 

2.773 

31 

348 

17 

*442 

• 556 

300 

0 

832*0 

31 

.9 

17.2 

315 

350 

2.574 

32 

258 

V 

*045 

.528 

350 

0 

901*0 

31 

.0 

17,6 

391 

200 

2*445 

31 

4.47 

18 

• 293 

• 582 

400 

0 

978*0 

31 

.8 

16.4 

467 

100 

2*3q7 

31 

056 

18 

• 750 

*604 

450 

0 

1038*0 

32 

• 2 

lfe.0 

559 

000 

2.236 

28 

571 

20 

• 134 

.705 

500 

0 

1118*0 

35 

• 0 

14.9 

488 

700 

2.413 

26 

667 

18 

• 868 

.708 

450 

0 

1086*0 

37 

.5 

15.9 

415 

200 

2*5g5 

27 

174 

18 

• 868 

• 694 

400 

0 

1038 • 0 

36 

.8 

15,9 

343 

000 

2*800 

29 

851 

18 

*634 

• 624 

350 

0 

980*0 

33 

.5 

16.1 

275 

700 

3*063 

29 

851 

18 

,634 

.624 

300 

0 

919*0 

33 

.5 

16.1 

211 

750 

3*388 

30 

675 

18 

•182 

• 593 

260 

0 

847*0 

32 

*6 

16.5 

153 

000 

3*825 

30 

030 

18 

• 182 

• 605 

200 

0 

765*0 

33 

.3 

16.5 

102 

300 

4*547 

27 

933 

17 

• 964 

• 643 

150 

0 

682*0 

35 

*8 

16.7 

61 

500 

6.150 

26 

882 

16 

• 575 

• 617 

100 

0 

615*0 

37 

.2 

18.1 

26 

650 

10*660 

24 

510 

16 

•216 

• 662 

§0 

0 

533*0 

40 

*8 

18.5 

12 

025 

19.240 

22 

989 

15 

• 306 

• 666 

35 

0 

481*0 

43 

.5 

19.6 

12 

025 

19*240 

22 

989 

15 

• 306 

• 666 

35 

0 

481*0 

43 

.5 

19.6 

7 

095 

31*533 

21 

368 

16 

• 484 

•771 

i'5 

0 

473*0 

46 

.8 

18.2 
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CHAPTER III 


FEASIBILITY OF CROSSED-BEAM APPLICATIONS TO COMBUSTIVE FLOWS 


by 

J. B. Thomison, Jr. and J. B. Stephens 
NASA - Marshall Space Flight Center, Alabama 


SUMMARY 


Group velocities of local areas of combustion in a turbulent flame 
have been optically measured by employing the crossed-beam correlation 
technique. Identification of velocity and photo-chemical reactions 
through correlation of emission fluctuations is discussed. 


III-l. INTRODUCTION 


The feasibility of applying the cross-beam technique was first 
developed by Fisher and Krause [ III- 1 ] in a subsonic jet. The method 
was extended to measure atmospheric wind speeds by Stephens, Sandborn 
and Montgomery [III-2] and later to measure optical fluctuations in a 
glow discharge described in Chapter II by Stephens and Thomison. The 
glow discharge measurements studied fluctuations of emission processes 
rather than extinction processes. Earlier tests depended on fluctua- 
tions in the transmission of an independent light source produced by 
changes in the optical density of a gas. These changes in optical 
density correspond to like changes in the mass density of the gas. 

The source for the wind tunnel measurements was artificial light, 
while that for the atmospheric tests was scattered sunlight. The 
tests described in this paper use the light emitted from a turbulent 
flame to study the combustion process of the flame itself. No outside 
light source is used. These tests demonstrate the feasibility of using 
the cross-beam technique to successfully retrieve velocity and spectro- 
scopic information from a combustion process. 
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Ill- 2. PHOTOGRAPHIC EVALUATION OF FLAME PROPAGATION 


Optical fluctuations in the body of a turbulent flame are due to 
turbulent eddies of burned and unburned gases (see Fig. III.l). These 
turbulent eddies may be defined as blobs of fluid which experience a 
coherent change of emissivity. In general, these blobs have a random 
velocity and after traveling a mixing length equal, on the average, to 
their diameter, they lose their separate identities [III-3], 

Before the actual cross-beam flame tests, the flame was photographed 
with a high-speed movie camera to determine the type of fluctuations pre- 
sent. The best visualizations of the flame motion were obtained at 
camera speeds in excess of 1000 frames per second. The movies were used 
primarily to determine the best position along the axis of the flame to 
observe with cross-beam photo-detectors. Qualitative velocity measure- 
ments are illustrated by the consecutive eight frames shown in Fig- 
ure III.l. The formation and motion of an emitting blob of fluid are 
observed at the right edge of the figure. The blob is moving to the 
right with a speed of approximately 5.0 m/sec and its axial scale is 
approximately 60 mm. 


II I- 3. CROSS-BEAM EXPERIMENTS 


Two mutually perpendicular lines of sight: were set to select the 
local point of investigation in the flame. All zero separation measure- 
ments refer to the selected position 10.1 cm downstream from the nozzle 
exit. To achieve separation measurements, the horizontal photo-multi- 
plier beam remained fixed, while the vertical beam was moved downstream 
as in Figure III. 2. The focal length of the optics was 11.3 cm with a 
beam diameter at an intersection of approximately 2 mm. The measure- 
ments were made primarily at 5250 A. 

A flame was achieved by feeding C 2 H 2 and 0 2 to the nozzle, with an 
oxygen-to-fuel ratio of 0.92. This ratio was set by adjusting the pres- 
sure controls of the feeder lines to 3 psig for C 2 H 2 and 1 psig for 0 2 . 
The actual oxygen and fuel mass flow ratio, calculated from the pressures 
assuming sonic speed at the nozzle exit, were 67.2 g/sec and 62.7 g/sec, 
respectively. The oxygen-to-fuel ratio was insufficient for complete 
combustion, and most of the burning occurred at the flame boundaries 
where ambient air was entrained. 

The actual experiments were made by feeding the outputs of the 
photo-multipliers directly into the PAR analog correlator to obtain a 
covariance reading. This on-line analog computer assumes that no D.C. 
shifts, trends, or gain instabilities were introduced by the electronics. 
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Further, the covariance thus obtained will become uncertain as soon as 
the correlation coefficient (correlation/mean signal power) drops below 
10 percent. 

Within the above assumptions, the analog computations of the PAR 
correlator should form an estimate of the temporal covariance function, 
which is defined as follows. 


C 

xy 


( T > ^ j A) 


...) = < 


1 1 > 
X y 


1 

T 


i (t;A) i (t- t; A) dt, 

x y 


( 1 ) 


where 


i 


x 


(t) 


yc-o 

T 

t 

T 


- instantaneous fluctuations of intensity from 
the undelayed channel 

- instantaneous fluctuations from the delayed channel 

- delay time set during data reduction 

- real time 

- total integration time. 


This covariance should be related to the statistics of local emission 
fluctuations. Let £ = q or £ = £ denote a point on one of the two 
beams. The emission fluctuations in the beam element AJ5 around this 
point are then described by the emission coefficient fluctuations 


K' U; t; A) 


A;t) 

A0 I(i;A;t) ’ 


( 2 ) 


where 

I(j}; ...) - radiative power flux through beam cross section at £ 

AE' “ emission fluctuation in beam element Af. , 

e 
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Crossed-beam theory [III-4] now indicates that the normalization of 

the A. C. -coupled signals with the D.C. mean values should provide a 

covariance estimate, which is equal to the point-area product mean value 

of the emission coefficient fluctuations K' around the beam intersection, 
— * 
x. 


G(|; t; A)— > 
2 \ 


< i i > 

x y 

i i 

x y 


= <K' (x,t;A) J J K' (x+6; t-x;A) dqd^ > , 

eddy 

correction 


(3) 


where 


= average D.C. intensity from beam x 
I = average D.C. intensity from beam y 

| = axial separation distance 

— ^ 

g = (|, T), (;) point coordinates in beam fixed frame. 


The last equation demonstrates that local emission variations in a 
flame can be retrieved by crossed-beam correlative techniques. The 
following pilot experiments were conducted to experimentally check this 
theoretical claim. 


III-4. VELOCITY MEASUREMENTS AND 
DETERMINATION OF COMBUSTION EFFECTS 


To obtain a velocity measurement in the flame, the two photo-multi- 
plier beams were separated along the flame axis, and the peak cross- 
correlation of these fluctuations revealed the average transit time for 
a turbulent blob. This transit time yielded the average velocity of the 
blobs over that distance. The results are illustrated in Figure III. 3. 
They show a decreasing velocity with an increase in the downstream beam 
separation. The speed of the eddy that was indicated by the noise falls 
into the speed interval shown in Figure III. 3. All photo-multiplier 
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records taken from single photo-multiplier time histories show a broad 
band power spectrum (see Chapter II). The signals could easily be 
resolved out of the photo-multiplier and flow noise, as apparent from 
the large values of the correlation coefficient for a zero beam separa- 
tion. The values were at least 20 percent, while some even exceeded 
50 percent (see table III.l), Further, the correlation curves were 
consistent and stable with repeating runs. When the separation was 
increased to 3.8 cm, the value of the correlation coefficient dropped 
below 10 percent, and no usable measurements were obtained. 

Additional tests were made to separate photo-chemical emission from 
heat radiation by studying the change of the cross-correlation as a 
function of optical wave length [III-4]. All these measurements were 
taken with intersecting beams. 

The results of this experiment showed distinct correlation peaks 
when changing the optical wave length. This would not have been pro- 
duced by heat radiation. Some of these emission peaks are centered 
close to optical wave lengths that correspond to known emission bounds 
for a C 2 H 2 - 0 2 flame [III-5], These results are illustrated in Fig- 
ure III. 4. This coincidence of the peak positions in crossed-beam 
spectra and published mean value emission spectra gives additional 
confirmation to the claim that local emission variations have been 
retrieved . 


III-5 . CONCLUSIONS AND RECOMMENDATIONS 


This research has demonstrated the feasibility of the use of the 
cross-beam correlation technique to successfully retrieve both velocity 
and spectroscopic information from a combustion. The next logical step 
would be to attempt similar studies in an actual rocket plume. 
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FIGURE III. 2. EXPERIMENTAL SETUP 
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FIGURE III. 4. PHOTO- CHEMICAL EMISSION IN A C 2 H 2 - 0 2 FLAME 



TABLE III.l. 


DATA FOR OBSERVED PHOTO- CHEMICAL EMISSION 


IN A C 2 H 2 - 0 2 FLAME 


X 


C xy { 0 ) 


C xx ( 0 ) C yy ( 0 ) 


6 ( X ) 


i -- 


0 


WAVE LENGTH 
ANGSTROMS 

2850.0 

2950.0 

3050.0 

3150.0 

3250.0 

3350.0 

3450.0 

3450.0 

3550.0 

3650.0 
3 750.0 

3850.0 

3850.0 

3950.0 

4050.0 

41 50.0 

4250.0 

4350.0 

4450.0 

4450.0 

4450.0 

4650.0 

4750.0 

4850.0 

4850.0 

4950.0 

5050.0 

5150.0 

5250.0 

5250.0 
5 350.0 

5450.0 

5500.0 


CORRELATION 
COEFFICI ENT 

0.238 
0.257 
0.261 
0.262 
0.246 
0.280 
0.335 
0.315 
0.303 
0.191 
0.201 
0 .331 
0.344 
0.302 
0.338 
0.322 
0.350 
0.334 
0.350 
0.381 
0.402 
0.442 
0.408 
0.396 
0.387 
0 .366 
0 .341 
0.391 
0 .492 
0 .506 
0 .331 
0 .223 
0.156 


EMISSION 

COEFFICIENT 

0.0163 
0 . 01 48 
0 .01 19 
0.0144 
0.0130 
0.0082 
0.0112 
0 . 0106 
0.0089 
0.0033 
0 .0041 
0.0065 
0.0049 
0.0054 
0.0052 
0.0047 
0 .0042 
0 . 0044 
0.0047 
0 .0060 
0.0060 
0.0075 
0.0075 
0.0085 
0 .0072 
0 . 0057 
0.0048 
0.0055 
0.0104 
0 .0080 
0 . 0103 
0. 0401 
0.06 70 
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CHAPTER IV 


TURBULENCE MEASUREMENTS BY THE CROSSED- BEAM TECHNIQUE 


by 

M. Y. Su 

Northrop - Huntsville, Alabama 
and 

F. R. Krause 

NASA - Marshall Space Flight Center 


SUMMARY 


An analytical relationship which approximates the probability 
density of wind component variations through the cross -correlation of 
the passive photometer outputs of the cross-beam detection system is 
presented. Specifically, the probability density is shown to be pro- 
portional to the cross-correlation of the first time derivatives of the 
time histories of fluctuating radiation intensities which are monitored 
by two non-parallel photodetectors. The direction of the measured wind 
component is normal to both lines-of-sight of the photodetectors. 

Special data conditioning and handling methods for evaluating the cross- 
correlation of the signal time derivatives are described, and their 
effects on the cross -correlation and corresponding spectral densities 
are discussed. The new analytical relationship and data handling techni- 
ques were applied to actual atmospheric crossed-beam field test data. 

The reduced optically and remotely detected probability density of wind 
component variations was found to be in reasonable agreement with that 
obtained by the conventional anemometer readings. Turbulent intensity 
can be easily evaluated once the probability density is measured. The 
optical technique can also be applied to other high- temperature, 
chemically reactive or supersonic turbulent flows such as rocket exhausts, 
jet plumes, and plasma in nuclear engines. The flow fields associated 
with these problems still defy accurate measurements by other presently 
available instrumentation. 
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NOMENCLATURE 


Symbol 

B 

d(i) 


E(£) 

f 


G AB^’ f) 




H(f ) 

h(T) 

H,( f ) 

h / t ) 


L„ 


m 


(t) 

(t) 


Definition 

bandwidth 

factor representing the degree of decay with respect 
to | 

function defined by equation (29) 
frequency 

cutoff frequency in a low-pass filter 
lower frequency in a time history 
terminate frequency in a low-pass filter 
upper frequency in a time history 
gain of S AB (|,f) 

gain of S^>(|,f) 

frequency complex frequency response 
filter transfer function 

frequency response of the Martin-Graham type roll-off 
low-pass filter 

transfer function of the Martin-Graham type roll-off 
low-pass filter 

fluctuating component of photometer A's output 
fluctuating component of photometer B's output 
integral time scale of atmospheric turbulence 
total number of digital samples 
total number of digital samples in m t ^ 1 class 
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Symbol 

P(i,T m ) 

P(U) 

R AB ^ T ) 

T p 


S AB^’ f ^ 



(l,f) 


T 

AT 


t 

At 

U 

U 

U 

m 

AU 

V 

X(f) 


x(t) 
Y (f ) 

y m^ 

y(t) 


Definition 

probability density of transit times in mth class 
probability density of wind component U 
cross-correlation of i^ and ig 

cross -correlation of samples from mth class 

cross-correlation of Si^/St and ^ig/^t 

spectral density of (| , t) 

^ 2 ) 

spectral density of (£,t) 

AB 

total length of photometer output 

length of each piece of photometer output 

time 

sampling interval of digital time history 
wind component normal to both light beams 
mean wind component 

wind component associated with mth class 
resolution in the wind component U 
wind speed 

Fourier transfer of x(t) 

kth sample of mth class of x(t) 

common signal in i^(t) 

Fourier transform of y(t) 

Hth sample of mth class of y(t) 
common signal in ig(t) 
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S ymb o 1 


Definition 


Greek Letters 


e(f,N d ) 


i 

S* 


a 


u 


T 


measure of error in the frequency response for using 
N d samples 

angle between the wind vector and the normal to both 
light beams 

normal separation distance between two beams 

normal separation distance between A and the imaginary 
beam B 

root mean square of the wind component U 
time lag between i^ and ig 


Operators 


< >i 


iAT 


AT 


( ) dt piecewise average of ( ) 


(i-l)AT 


( ) accumulative average over m piecewise averages 

m 


( ) 


m 





confidence level of 
accumulation average 


&R 


m 


+Jm 


1 


a a 
xm ym 



R (j) 
m VJ / 



relative statistical 
error of correlation 
or covariance 
functions . 
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IV- 1 . INTRODUCTION 


Remote optical detection techniques are being developed by the 
Marshall Space Flight Center to provide a means for continuous monitor- 
ing of atmospheric wind and turbulence. Successful detection of most 
probable wind speed at an altitude of about 100 meters by these techni- 
ques has been reported [IV-1 5 -2], 

The capability of these techniques has been extended to atmos- 
pheric turbulence measurements by optically retrieving the probability 
density of wind component variations. Direction of the wind component 
is normal to both lines-of-sight of the photometers employed. The mean, 
rms , and other higher moments of wind component variations may be easily 
evaluated from the probability density. Preliminary results on this 
progress with emphasis on the signal/noise analysis and electro-optical 
system were reported in reference IV-3. 

The optical technique essentially consists of two photodetectors to 
monitor the optical fluctuations of natural scattered radiations which 
result from both spatial and temporal variations of atmospheric inhomo- 
geneities. Each photodetector output will thus contain the information 
of the atmospheric inhomogeneities along its entire field of view. The 
cross-correlation of these two photometer outputs should reveal the 
characteristics of the local optical disturbance and atmospheric motion 
which are common to both fields of view of the two photodetectors. To 
have an optimum spatial resolution the two lines-of-sight of detectors 
should be so oriented that they are normal (but separated) to each other. 
This is the reason for calling this optical technique the crossed-beam 
technique. Because of the very low signal-to-noise ratio associated 
with the remotely detected radiation fluctuations, special data condi- 
tioning and processing and statistical analyses are needed for evaluating 
the cross-correlation of the photometer outputs and their time derivatives. 

This paper presents the theoretical basis for optical approximations 
of turbulence levels, together with the necessary additional considera- 
tions on data conditioning and processing. In section IV-2 , a theoret- 
ical relationship between the cross -correlation of the first-time 
derivatives of the photometer outputs and the probability density of 
the wind component variations is derived. This wind component is normal 
to both lines-of-sight of the photometers. Use of the time derivative 
of the photometer outputs instead of the outputs themselves is one of 
the essential features which differ from the previous crossed-beam 
detection of the most probable wind component. 

Section IV-3 describes the special data conditioning and processing, 
and discusses their effect on the cross -correlations and corresponding 
spectral densities. These considerations are very important for dealing 
with cross -correlating time derivatives of random time histories. 
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Section IV-4 presents some experimental and computed results from 
applying the theoretical relationship and data handling methods to atmos- 
pheric crossed-beam field test data. The optically detected probability 
densities of the wind component variations were then compared with those 
independently measured by the conventional anemometers mounted on a 
close-by meteorological tower. The minimum separation of the crossed- 
beam arrangement was set directly upstream of the anemometer about 
60 meters above the ground. 

The authors wish to thank the following persons for their contribu- 
tions to this paper: Professor V. Sandborn, Colorado State University, 

for carrying out the atmospheric crossed-beam field tests; Messrs. J. 
Jones and R. Graham, Computation Laboratory, MSFC, for digitization of 
test tapes, cross-correlation computation and the program for digital 
filtering; Mr. J. Pooley, Northrop Corporation, for programming an 
updated version of the error accumulation and spectral decomposition 
subroutines; and Dr. W. H. Heybey, MSFC, and Dr. L. N. Wilson, IITRI, 
for their very helpful and constructive critical reviews and comments 
on the extension of the crossed-beam concept to turbulence measurements. 


IV- 2. PROBABILITY DENSITY OF 
WIND COMPONENT VARIATIONS FROM OPTICAL CORRELATION 


Wind speeds and directions vary, in general, fairly rapidly. To 
specify completely the statistics of such wind characteristics, either 
an infinite series of correlation functions of all orders of wind fluc- 
tuations, or an infinite series of all orders of probability densities 
of wind fluctuations will be needed. However, in usual meteorological 
applications, only limited quantities are measured, i.e., the first few 
moments of the wind fluctuations such as the mean wind speed, and the 
rms wind speed fluctuations. 

When the wind direction is constant but unknown, and the speed is 
varying, it has been shown [IV-4, -5] that both the wind direction and 
speed can be simultaneously measured by the crossed-beam system. When 
the wind direction varies, only the wind component normal to the two 
lines-of-sight of detectors can be measured. The reason is that the 
crossed-beam measurement gives the transit time of the optical disturb- 
ance between two beams only, while the location of the transit path and 
its length between two beams cannot be determined unless the wind direc- 
tion remains constant over the period of measurement (see Figure IV. 1). 
Fortunately, the fixed position of two light beams gives a uniquely 
determined minimum separation between two beams which is in the direction 
normal to both beams. Thus, the combination of the measured transit time 
and the shortest beam separation always determine uniquely the wind com- 
ponent normal to both beams. Since, in most meteorological conditions, 
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the wind direction does vary in the minimum required period of crossed- 
beam measurements, we shall deal with this wind component throughout this 
paper . 

The first-order probability density of the wind component can be 
measured by the crossed-beam system to a high degree of approximation. 
This approximation results from some relatively weak assumptions intro- 
duced in order to establish an explicit analytical relationship between 
the cross-correlation of optical fluctuations and the probability of wind 
component variations. Also, because of the variation of the wind direc- 
tion, the resolution of the altitude of the prevailing wind will be 
questionable from one pair of crossed beams. However, if more than two 
pairs of crossed beams are used, which are so oriented that mutually 
perpendicular wind components are to be measured, then the wind direction 
can be approximately determined also. 

The basic geometrical arrangement of the crossed-beam system is 
shown in Figure IV. 1. The various fields of view of two ground-based 
photometers, A and B, are oriented so that their lines -of-s ight differ 
by 180 degrees in the azimuth angles. Such an orientation will give the 
wind component both horizontal to the ground and normal to both beams. 

The wind component will be the one used for this discussion, unless 
otherwise specified. 

Atmospheric inhomogeneities, such as the variations in density of 
aerosol, various gases and water vapor, will cause both spatial and 
temporal variations in radiation processes. These radiations include 
light scattering, absorption, and emission. Thus, the characteristics 
of local light fluctuations should reflect the characteristics of the 
atmospheric inhomogeneity (or turbulence). The atmospheric inhomogene- 
ities are transported according to the general rules of fluid dynamics. 
The lifetime of the atmospheric inhomogeneity, which is a measure of 
the average time lapse for the inhomogeneity to retain its character 
(or signature), is finite; so is the lifetime of the optical fluctuation, 
if the light source can be considered to be relatively constant. The 
optical detectors A and B thus record the time histories of accumulative 
radiation outputs within their entire field of view, respectively. Fluc- 
tuations of the photometer output, iy^(t) and ig(t), therefore contain 
information on the atmospheric motions along the entire line-of-sight. 

The mean value of the photometer outputs is not of concern here. The 
fluctuating part is separated from the total recording signal by A.C.- 
coupling electronic circuitry. 

When the two beams actually intersect each other at a certain height 
(Figure IV. 2), the local optical fluctuations radiating from the inter- 
section region will result in some common signals in the outputs of the 
detectors. The common signals may not be exactly the same, but they are 
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related, so that the photometer outputs may be expressed as 


i A (t) = x(t) + i Ajn (t) 
i B (t) = y(t) + i B>n (t). 


( 1 ) 


Here, x(t) and y(t) are the common signals mentioned above, while i^ n (t) 
and iB,n(t) are light fluctuations outside the intersection region, ’in 
contrast to desired common signals, they will be called noises due to 
transmission, or flow noises. By the nature of the atmospheric inhomo- 
geneities, all the terms in equation (1) are random time histories. The 
most powerful method available to retrieve such common signals is the 
correlation technique, which is commonly used in the communication theory 
and random vibration analysis. A correlation function is computed by 
adding the products of two time histories, i^ and i B . That is, 


e ab« 


<i A ( t ) i B ( 1:+ t) > s lim Y 

T— >00 


l A (t) 


i B (t+T) 


dt . 


( 2 ) 


Substituting equation (1) into (2) yields 


R AB (t) = <x(t)y(t+T)> + <x(t) i 0 n (t+ T )> + <i A n (t)y(t+-r)> 


+ <i (t) i (t+ T )> . 
A,n B.n 


(3) 


The first tertr^in the right-hand side of the above equation is the cross- 
correlation between two common signals. The other three terms are the 
cross-correlation among statistically unrelated light fluctuations; thus, 
they will vanish if the integration time T is long enough. That is, the 


time integration will suppress noise 
tion, the integration time T is always 
due to this finite T will be discussed 
ficiently long integration time, which 
histories, equation (3) can be written 

r ab ('0 ~ < x (t)y (t+x)>. 


components. In actual applica- 
finite. The uncertainty (or error) 
further in section V. For a suf- 
is also the length of the time 
approximately as 


(4) 
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This means that as long as one deals with such product mean values, one 
can disregard all light modulations that are not common to both fields 
of view of the two detectors [IV-6j. 

To measure the wind component the two beams have to be separated. 
The shortest separation distance is denoted by |. The common signals 
in the two detectors now are caused by the convected atmospheric inhomo- 
geneities which happen to traverse both light beams. The time lapse 
which an inhomogeneity takes to travel from one beam, say A, to the 
other beam, B, will be called the "transit time." The distance between 
the two beams of the traversing path CD will be called the "transit 
distance." Associated with the transit path is the transit height, 
which is the altitude of the transit path above the ground (see Fig- 
ure IV. 1 ). When the wind speed, V, and direction, 9, are fixed, the 
associated transit time can be obtained from the cross-correlation 
curve (equation ( 4 )). It is the time lag associated with the maximum 
value (or peak) of Rab( t )> denoted by T pe ak> 0n t * ie other hand, the 
transit distance cannot be determined unless the wind direction is 
known. Without knowing the wind direction, the best one can deduce is 
the component of the transit distance in the direction normal to both 
light beams, which is simply the shortest distance Thus, the 
relationship between the wind component U and the time lag Tp ea k from 
the cross ed-beam measurement is 


CD £/cos 9 

U = V cos 0 = cos 0 = — cos 0 

peak T peak 

or 

U = ■ 7 L - . ( 5 ) 

peak 

Here, the distance (• is fixed by the crossed-beam arrangement, which is 
completely known. Hence, equation ( 5 ) allows us to measure the wind 
component through the remote detection of the optical fluctuations and 
the cross -correlation of these fluctuations if the wind speed and direc- 
tion remain constant. 

The more realistic meteorological situations are those where the 
wind speed and direction both change with time. Under such situations, 
the maximum information that can be retrieved from the cross-correla- 
tion of the optical fluctuations is the first-order probability density 
of the wind component variations . 
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We shall now derive the relationship of the cross -correlation of 
crossed-beam signals to the probability density of the wind component 
variations. As discussed above, all the noise components in the two 
photometer outputs i^(t) and ig(t) will be suppressed through cross- 
correlation. We need only to consider the two time histories of com- 
mon (desired) signals, x(t) and y(t), without loss of generality. Let 
the time histories of length T be digitized with a sampling interval, 
At. The total number of digital samples for either time history will 
be denoted by N: 


N 


X 

At * 


These digital samples and digital time lags will be denoted as follows: 


= x (t = n A t) 
y n = y(t = n&t) 


and 


( 6 ) 


T 


l 


SL At, 


where n = 1,2, ..., N and H = 0,1,2, L«N. 

Since the wind component may vary over the period of measurement, 

T, the associated transit time may also vary correspondingly. Without 
loss of generality, we shall assume that the beam A lies in the upstream 
of beam B. Each digital sample x n will then have a corresponding digital 
sample y n +£. They are attributed to the same atmospheric inhomogeneity 
which takes the time lag = £ At to travel from some point in beam A 
to another point in beam B. The transit path between these two points 
is not necessarily straight. It can be of any shape. The corresponding 
wind component normal to both beams is |/t^. The cross -correlation 
(equation (4)) can now be written in the summation form as 

N 

R AB ( V = N X X n Y n+i* (7) 

n=l 


It should be noted that the value of the summation is invariant under 
any different ways of grouping the terms x n y n+ ^ in the summation. For 
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the purpose of retrieving variations of the wind component, the digital 
samples x n may be classified with respect to their associated transit 
times. Let the class interval for transit times be the same as the 
sampling interval At. The mth class of transit times, denoted by t 
will be defined as 


t : m At ^ t < fm+1) At, (8) 

m m 


where m = 0,1,2, . . . , M. M is any integer that is large enouth to cover 
all possible transit times of interest. The corresponding mth class of 
wind components denoted by U m will be 


U : 
m 


I 

(m+1) A t 


< U 


m 


m A t 


(9) 


Now group those digital samples x n according to the classification of 
equation (8). Further, denote the kth sample in the mth class in the 
ascending order of time by 


x m = x(t = m k At) (10) 

k 

and the total number of digital samples belonging to this class by N m . 
Thus , 

M 

N = ) N . (11) 

A_j m 
m=0 

The probability density of transit times of the mth class over the period 
T of measurement can be expressed as 


P(£,T m ) A t = PU, 


m A t g t < (m+1) 
m 


At] At 


N„ 


( 12 ) 
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The corresponding probability density of the wind components P(U m ) can 
be obtained from the above probability density by a transformation of 
variables from the transit time to the wind component. 


p <v 


= p 


< U s — = — 

(m+1 ) A t m m A t 


P (I,T ) 


i dT mi 


m/ >dU, 


p a.%) 


(13) 


m 


With the above classification of digital samples and definitions of the 
probability density, the cross-correlation of the optical fluctuations 
from the photometers can be expressed as 

N 

R AB ( V = U x - y n +i 

n=0 


M ^m 



m=0 k=l 


Xm k Ym k + -« 



m=0 k=l 


(14) 


Let us define the term in the parentheses by 

^m 

R » <6 ’ T ^ 

k=l 


(15) 


By substituting equations (12) and (15), equation (14) can be written as 


M 


R AB (l ’ T i ) 


N, 


m 
m=0 


p (^ T m ) KmU.Tg) At 


(with At = At) 
m 


(16) 
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Clearly, the function R m (£,xj;) represents the cross -correlation of the 
mth class of digital samples. Equation (16) simply states that the total 
cross-correlation R^gC^x^) is equal to the sum of M + 1 cross-correla- 
tion of each individual class of samples with the associated probability 
density of transit times during measurement T as a weighting factor. 

In almost all actual applications, both functions P(^,x m ) and 
Rm(£> T m) can be assumed to be continuous with respect to x m ; thus, the 
above summation can be written in the integral form as the increment 
Ax m becomes very small (i.e., Ax m 0) . 


R AB^’ T p 


P(|,t ) R (|, x .) d x . 
m m i m 


(17) 


Strictly speaking, this limiting process Ax m -> 0 cannot be applied, 
since this would require the knowledge of the infinitesimal wind com- 
ponent fluctuations. It is known that the highest wind fluctuation 
measurable by the crossed-beam technique is dictated by the beam 
separation distance | and the maximum prevailing wind component. 
Specifically, the highest resolution in the wind component is given 
by 

AU = , (18) 

T max 


where x max is the maximum time lag used in R^g. This is the lowest 
possible wind component which can be measured by the crossed-beam 
system, i.e., 


min(U) 


= AU = 


_i_ . 

T 

max 


(19) 


The highest frequency fluctuation of the wind component is given by 

max (f ) = , (20) 

u / rt | 
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where U max is the maximum possible wind component in that particular 
measurement under consideration. However, we applied the limiting 
process in the sense of numerical equivalence between equations (16) 
and (17), which offers a high degree of approximation. Equation (17) 
is an integral equation for the desired (unknown) probability density 
of transit times. The simplest possible way to solve for P(|,T m ) is 
to uncouple interactions between P(£,x m ) and R m (|,x^) so that the 
probability density can be expressed in terms of the experimentally 
accessible cross-correlation. The desired uncoupling is feasible with 
suitable data conditioning on the photometer outputs under some rela- 
tively weak assumptions regarding atmospheric turbulence and wind 
characteristics. This will now be discussed. 

We shall first make an assumption (I) that samples in each class 
will decay by the same amount in their turbulent characteristics over 
the same beam separation and that the rate of decay may differ in dif- 
ferent classes of samples. According to Bachelor [IV-7], the decay of 
relative turbulent energy, which is defined as the ratio of the mean 
square of speed fluctuations to the square of the mean speed, is linearly 
proportional to the transit distance. The turbulent fluctuation in the 
photometer outputs is a mixed outcome of fluctuations both in the flow 
speed and the extinction coefficient of the flow medium. Nevertheless, 
the assumption would be at least the first order approximation for a 
transit distance of the order of several integral length scales of 
atmospheric inhomogeneities. The second assumption (II) is that the 
number of samples in each class, which is classified at the upstream 
beam A, is the same if the classification were made at the downstream 
beam B. This assumption is really a restriction which excludes the 
consideration of wind component fluctuations with frequency higher than 
max [f u ] (see equation (20)). By using these two assumptions, the cross- 
correlation of mth class of samples (equation (18)) can be rewritten as 

V&.Tg) = D(£) R(i* = 6 - u m x r o) (21) 

or 

VI, *e) = R (£* = K 1 - -“)> 0). (22) 

l m 


where D(|) is a factor representing decay of turbulence with respect to 
the separation distance, while the space-time cross-correlation R m (i,x^) 
is replaced by an equivalent spatial cross-correlation R(|*,0) with 
imaginary translation of the downstream beam B closer to the upstream 
beam A by the amount of U m x for compensating the convective motion of 
the atmospheric inhomogeneity (see Figure IV. 3). 
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Next, substituting the approximate relationship (equation (22)) into 
equation (17) yields 


R ab U.^) - D(t) 


P(S.%) R(6* - 6(1- yb, o) d-T, 

l m 


(23) 


What has been accomplished in the above equation compared with the equa- 
tion (17) is that the dependence of the functional form of R ra has been 
eliminated from the integrand. However, P(|,T m ) and R(|“,0) is still 
coupled. It should be noted that the spatial cross-correlation R(|'",0) 
reaches its maximum value for |* = 0, i.e., when the imaginary beam B* 
intersects beam A. A certain correlation will persist for 0 < , but 

will be smaller than the scale of the atmospheric inhomogeneities common 
to both light beams. In other words, a certain correlation exists if 

I T m T il = I hi/ 1 ' T''>\ 

is less than the time lapse for the atmospheric inhomogeneity traversing 
its own length. This time lapse is usually defined by 


L 


T 


1 

R(£* = 0 , o) 


oo 

J R(i*,0) 

o 


(24) 


which is normally called the integral time scale of the atmospheric 
inhomogeneity. For the time being, we shall introduce a third assump- 
tion (III) that the integral time scale is much smaller than the half- 
width of the probability density of the transit time fluctuations (see 
Figure IV. 4). That is, 


L t « T 


rms 


= 


(25) 


where t 2 and i 2 are the two transit times such that 
P(£,Ti) = P(6,t 2 ) = e“ l/s max (P (| , t) ) . 


(26) 
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The justification of this assumption will be discussed after equation 
(36). But let us examine first how this assumption will affect equation 
(23). This assumption essentially implies that the spatial cross-corre- 
lation R(|*) behaves like a Dirac delta function. The equation (23) can 
be approximately written as 


R AB^’ 


D(|) PCl.Tg) 


R(e = 6 ( 1 - 


■~)>°) d ^m- 
l m 


(27) 


Changing the variable in the integration from x m to yields 


dx 


r AB ( ^ t P = D ^) 


PCi,^) 


R(f ) 


d| 


d|' 


= D(|) P(6, Ti ) 


MO. 


(i- 


10 


d|* 


= E(|) T(g Pd.Tg), 


(28) 


where 


e(i) 


Mil 


00 


— 00 


R(f ') 

(l- £-) z 


a. 

dr. 


(29) 


Thus, the probability density of transit times can be written as 


Pa.Tg) 


r ab^ , t P 

E(£) T i 


(30) 


The dependence on the factor E(|) can be further eliminated by normaliza- 
tion, since it is independent of the transit time. That is, 
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P(i,Tg) 


( 31 ) 


max(P(|, tJ) 


R AB^’ T P / ' T l 

max(R AB (|, t )/t^) 


The normalized probability density of transit times can be approximated 
by the measurable cross -correlation of two photometer outputs. Physi- 
cally, it will be more useful to have the probability density of wind 
component variations. This can be accomplished by transforming the 
variable from the transit time % £ to the wind component U through the 
following relationship: 


U 



(32) 


The normalized probability density of wind components (normal to both 
light beams) P(U) can be approximated by using equation (30). 


P(U 


■f) = ) 

T i * 


_ h 

1 dU * £ 


P(&,Tg) 


R AB ( ^’V 

£ E(|) 


(33) 


Again, by applying the normalization procedure to the above expression, 
we obtain 


p(u= M 

max (P (U) ) 


r ab ( £'V 


max[T i r ab ( 6 *V ] ‘ 


(34) 


Alternatively, one may apply the characteristics of the probability 
density, i. e. , 


00 

y'p^.Tg) = i > ( 35 > 

-00 
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to equation (30), and then substitute in equation (33): 


P(U= — ) = 


h r ab ( s - t P 


I / — r r > „(p.t' ) d T ' 

J r SL 


( 36 ) 


This implies that the probability density of wind components normal to 
both light beams can be approximated by the measurable cross-correlation 
of photometer outputs. The above optical approximation is valid under 
the three assumptions made which are summarized as follows: 

(I) Samples of common optical fluctuations in each class with 
the same transit time will decay by the same amount in 
their turbulent characteristics over the same beam 

s eparation. 

(II) The number of samples in each class, which is classified 
at the upstream beam, is the same if the classification 
were made at the downstream beam. 

(III) The temporal integral scale of atmospheric inhomogeneities 
is much shorter than the half-width of the probability 
density of the transit times. 

The first two assumptions are fairly realistic approximations for atmos- 
pheric turbulence, provided the transit distance does not exceed about 
200 meters or the transit time does not exceed about 20 seconds. On the 
other hand, the third assumption is usually not met in the atmosphere. 

The integral length scale of atmospheric inhomogeneity is between 10 to 
100 meters long in the atmosphere [IV-8]. However, a shorter optical 
disturbance may be expected if we cross -correlate the first time deriva- 
tives c3x(t)/c)t and By(t)/cSt of the common signals x(t) and y(t) instead 
of the common signals themselves. Since we shall employ the cross -corre- 
lation technique, these two time derivatives can be simply replaced by 
the time derivatives of the photometer outputs i)i^(t)/c)t and <3ig(t)/c)t. 
Physically, the integral length scale of the time derivatives is the 
so-called microscale of optical disturbance [IV-9]. This microscale is 
usually at least one order shorter than the average integral scale. 
Mathematically, the time derivative is the time rate of change of optical 
disturbance; thus, significant values of the derivative may be expected 
only when a large change in the optical disturbance occurs. Heuris tically , 
if one may simply visualize an atmospheric turbulent eddy as an elongated 
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body with its major axis along the flow direction, then the time deriva- 
tive will be directly proportional to the slope of the eddy. Therefore, 
only when the leading or trailing edge of the eddy passes through the 
light beam may a significant optical fluctuation occur. Therefore, the 
unrealistic assumption (III) can now be replaced by the following real- 
istic assumption (III-A), if the time derivatives of the photometer out- 
puts are used to evaluate the cross-correlation. 

(III-A) The temporal microscale of atmospheric inhomogeneities 
is much smaller than the half-width of the probability 
density of the transit times. 

With this consideration, equations (34) and (36) should then assume 
the following forms: 


P(U= ^) 

Z 

max (P (U= -^-) ) 

T i 


R AB >(t 'V 

( 2 ) 

maxf^ R ab U,t £ ) 


(37) 


and 


where 


P(U= f-) = 

z 


T R^ 2 -* ( £ T ) 

Z AB U, V 






1 

T 


Si^(t) Sigdjt+x) 


St 


St 


dt. 


(38) 


(39) 


For statistical stationary processes, the above expression is equivalent 
to the negative of the second derivative of the cross-correlation; i.e.. 


r ab ) 


S 2 

St) 


R AB^ ,T i^ 


(40) 
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Once the probability density of the wind component variations is obtained 
from the cross-correlation of the time derivatives of photometer signals, 
all the central moments of the wind component variations can be easily 
calculated. The two most important moments are the mean wind component, 


U = 


00 



U P(U) du, 


and the mean square value (or variance) of the wind component variations, 


a 


2 

u 


00 

f (U-U ) 2 

— CO 


P(U) du. 


The ratio of the root-mean-square value to the mean wind component 
(a /U) is a common measure of atmospheric turbulent intensity. 


IV- 3. DATA CONDITIONING AND PROCESSING METHODS 


It was shown in the last section that the time derivatives of photom- 
eter outputs are needed for evaluating cross-correlation in order to 
retrieve remotely the probability density of wind component variations. 

For stationary atmospheric conditions, this will be equivalent to require 
evaluating the second derivative of the cross -correlation of the photom- 
eter outputs themselves. The atmospheric phenomena are known to be 
intrinsically nonstationary, but for qualitative understanding of the 
problems involved in the data conditioning and handling, we would assume 
for the moment that the time histories of photometer outputs are stationary. 
Under this assumption, equation (40) holds. Taking the Fourier integral 
transform of the equation yields 



( 1,0 


f2 s ab«' £ >' 


where 


,(2) 

’AB 


(i,f) 


00 



- CO 


(5 ’V 


i2itf 


T v 


(41) 


(42) 
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and 


s ab<£> £ > = / r ab«>V el2 ” £T< %• 


AB 


(43) 


Here, S^g(£,f) is the cross-spectral density of the photometer outputs, 

(g) 

while (|,f) is the cross-spectral density of the first time deriva- 

i\.ij 

tives of the photometer outputs. These two spectral densities are 
generally complex functions. For easier discussion, we would rather 
deal with their gains which are always real and defined by 




.( £ ) 

"AB 


AB 


and 


(44) 


G AB<^ £ > ' <45) 

where the asterisk indicates the complex conjugate. Hence, equation (41) 
can be rewritten as 


- £ 2 g ab ( ^ £ >- 


(46) 


This equation states that the gain from the time derivatives of the 
photometer outputs is equal to the gain from the outputs themselves 
multiplied by the frequency square. 

In the atmosphere, the gain of the photometer outputs G^g(^,f) 
follows approximately f -2 in the range of energy-bearing eddies. This 
would imply that the gain from the time derivatives would approximate a 
white noise spectrum. It is also well known that the wider the bandwidth 
of a white noise spectrum, the narrower the first zero crossing of the 
corresponding correlation function. Thus, to better approximate the 
assumption (III-A) , we would like to have the bandwidth of the photom- 
eter outputs as wide as possible. However, because of the following two 
factors, the bandwidth of the signals (i.e., the photometer outputs) 
should be suitably limited. First, the common signal-to-noise ratio of 
the photometer outputs with natural scattering radiations as light 
sources is very low. This ratio is usually of the order of 0.4 or less 
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[ IV- 3 ] . It follows that the correlation coefficient of the i^ and i B , 
i.e. , 


<X A 1 B > 

<i 2 > X / 2 <i2>l/2 
A £> 


(47) 


is of the order of (0.4) 2 = 0.16 or less. Such a small correlation coef- 
ficient could not be accurately evaluated by commercially available analog 
correlators. A digital computer is required to achieve a more complete 
noise cancellation by integrating over very long records without dynamic 
time placement errors. The f 2 weighting of the high frequency noise which 
is introduced by differentiating the time series also requires an addi- 
tional low-pass filter with a sharp roll-off at the high frequency limit 
of the common signals. Secondly, in using digital computers, the time 
derivative may be approximated by the finite difference of two consecutive 
digital time samples. The quantization errors in the individual samples 
resulting from the analog/digital conversion will further propagate to 
the time derivatives [IV-11], The main effects of the above two factors 
are significant only in the higher frequency portion of the cross-spectral 
density. 

To reduce the quantization errors in the digital time histories, the 
sample reduction method is employed. In essence, this method consists 
of two steps: (1) sampling the analog time history at an interval by 2m 

times smaller than that required by the Nyquist's criterion, i.e., 

At = l/2(2f u ), where f u is the upper frequency of the time history meas- 
ured in Hertz and m is any positive integer equal or greater than 1, and 
(2) averaging every m consecutive digital samples to form a new digital 
time history with an equally spaced time interval of 1 / (4f u ) . The theory 
[IV-11] shows that, on the average, about 30 percent reduction in the 
quantization error may result by reducing five samples to one on a time 
history having an f -2 spectral density. Reductions for other spectrum 
shapes are also given in reference IV-11. 

To determine optimum injection of f 2 weighted high frequency noise 
components, we have various standard frequency settings for the digital 
filters. Both center and roll-off frequencies were changed independently 
of each other. In applying the digital filtering, any desired filter 
gain and phase functions can easily be accomplished by using a suitable 
transfer function without worrying about practical limitations in 
electronics . 
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Specifically, let us consider an analog time history x(t) and let 
its Fourier transform be X(f), i.e., 


X(f) 



-i2jtf t 

e 


dt. 


“00 


(48) 


This function represents the frequency decomposition of the original 
time history. In a broad sense, the filter is meant to be an operator 
on X (f ) such that X(f) will be weighted according to a desired fashion 
in a certain range of frequency. Mathematically, this operation can be 
expressed as 


Y (f ) = H(f ) X (f ) . 


(49) 


Here, H(f) is the function characterized by the desired operation and is 
called the (complex) frequency response function of the filter, while 
Y(f) is the frequency decomposition of the output filtered random time 
history. The frequency response function can be broken into the real 
part H r (f) and imaginary part H(f), i.e., 


H(f) = H (f) + iH (f). (50) 

The absolute value | H(f ) | of H(f) is called the gain of the filter, 

| H(f ) | = [Hj(f) + H?(f) ] 1 / 2 , (51) 


while the arc-tangent of the ratio H^(f)/H r (f) is called the phase func- 
tion of the filter, 


Hi(f) 

0(f) = arc tan ( — ^ y) . 


(52) 
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Denote the inverse Fourier transforms of H(f) and Y(f) to the time 
domain by 


h(t) 


CO 



exp(i2rtft) dt 


( 53 ) 


and 


y(t) 


co 



exp(i2jtft) dt. 


-00 


(54) 


Taking the Fourier transformation of equation (49), and then sub- 
stituting the above expressions, we obtain 


y (t) 


h(t' ) x(t-t' ) dt' 


(55) 


This equation states that the filtering as defined by equation (49) 
is equivalent to an integral transform of the original time history x(t) 
with the weighting function h(t). This weight is called the transfer 
function of the filter. The output of this transformation y(t) is usually 
called the filtered time history. 

For an ideal low-pass filter with bandwidth B = f u , the frequency 
response function is simply 

r 1 > l f I = B 

H(f ) = H(-f) = < (56) 

U, I £ I > B. 

The corresponding transfer function will be 

h(t) - B (57) 

However, because of the Gibb's phenomenon, we cannot use the transfer 
function of equation (57), directly. This phenomenon essentially states 
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that oscillations exist in the approximating frequency response func- 
tion near the discontinuities of the ideal one when one approximates 
an ideal frequency response function having one or more jump discon- 
tinuities with a truncated Fourier series. 


To avoid this difficulty, the frequency response function is 
usually constructed as a continuous function. There are numerous ways 
to accomplish this. What we shall adopt below, one of the simplest 
types, is called the Martin-Graham type roll-off filter. Its frequency 
response function is defined for a low-pass filter [IV-12, -13] as 
foil ows ; 


r 


i. 


v« 



[1 + cos rt(f-f c )/Af ] , 
[1 + cos jt(f+f c )/Af] , 




0 , 


Id s f c 

f C < f < frp 
f T < f < - f c 

|f| § f T 


(58) 


where f c and f-j are the cutoff frequency and terminate frequency, 
respectively, while Af is the frequency interval between these two 
frequencies, i.e., Af = fx - f c . The inverse Fourier transformation of 
equation (58) yields the transfer function for the Martin-Graham low- 
pass filter: 


s in (2 irf t) -t sin( 27 tf T t) 

h (t) = — . (59) 

1 2«t[l - 4(Af) 2 t 2 ] 


The theory for filtering and the associated equations are formu- 
lated in terms of analog time histories. To apply them to discrete 
digital time histories, all integrations need to be rewritten in sum- 
mation forms. The frequency contents of the time history should also 
be known to set a suitable sampling rate for digitization following 
the Nyquist criterion. In actual application, the limits of integra- 
tion in equation (55) will be finite. 

Obviously one would like to limit the integration limits T as 
much as possible in order to minimize the required computer time. An 
explicit relation for the minimum integration time is now derived from 
the characteristics of both the time histories and the desired filter. 
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Let us denote the upper frequency in the time histories by f u , the desired 
cutoff frequency of the filter by f c , and terminate frequency by f>j. The 
upper frequency will determine the minimum sampling interval by Nyquist's 
or Shannon's criteria as 


At 


1 

2f ‘ 
u 


Then equation (55) can be written in summation form as 


N, 


y k = At z, Vk-f ’ 

£=-N, 


(k = N d , N d + 1, , (N - N d ), 


(60) 


where 

y k = y(k A t) 
h £ = h(£ At) 

~ x[(k-^) At] (61) 

and 


N = T/At. 


The limit of summation ±N d is dictated by the tolerance of error in the 
filter frequency response. Let us define the measure of error [IV-12] 
as 


e(f,N d ) = H(f ) - H(f,N d ), 


(62) 


where H(f,Nd) is the N d -term truncation of the infinite series represen- 
tation of H(f). It has been determined empirically that for the toler- 
ance 


0.57o < e < 1 % 


(63) 
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the minimum number required is given by 


5f 

mm ( N d ) = f ) > 

T c 


( 64 ) 


with 


This shows that the sharper the roll-off rate in the filter, the 
greater the number of digital samples needed for a fixed accuracy toler- 
ance. Figure IV. 5 shows the phase and gain function of the Martin-Graham 
low-pass filter with f c = 0.8 Hertz, f^ = 1.0 Hertz and = 39. The 
figure clearly shows that there is only negligibly small phase distor- 
tion and that the gain is almost exactly as desired. 


IV- 4. EXPERIMENTAL AND COMPUTED RESULTS 


To verify the theory that the probability density of wind component 
variations can be approximated by the cross ed-beam technique, and to test 
the applicability of the data conditioning and handling methods described 
in the last two chapters, some computer-processed results of actual 
crossed-beam test data recorded in the lower atmosphere are presented in 
this chapter. 

The field tests were carried out by the Colorado State University at 
Boulder, Colorado. A schematic crossed-beam test arrangement is illus- 
trated in Figure IV. 6. The detailed test site topography, geometry of 
the photodetectors, and meteorological conditions were reported in refer- 
ences IV- 1 and IV- 2. 

The light source is the natural scattered radiation in the visible 
spectral range. A block diagram of the crossed-beam detection system is 
shown in Figure IV. 7. Essentially, the incident radiation in the field 
of view of the telescope (0.5 degree) is focused on a photodetector 
(silicon diode). The radiation power is converted into voltage in pro- 
portion to its power. The fluctuating voltage is electronically separated 
into the d.c. and a.c. components with a time-constant of 100 seconds. In 
addition, a low-pass and high-pass filter of 0.01 Hz and 3 Hz, respectively, 
are set in the recording chain. Both a.c. and d.c. signals, together with 
the wind speed and direction obtained by anemometers, are recorded on an 
analog magnetic tape. More details of the optics and electronics, and 
associated signal/noise analyses were discussed in reference IV-3. 
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We shall present and discuss two test runs (I and II). The perti- 
nent information is summarized in Table IV- 1. The top view of the 
photometer arrangements and the probabilities of wind directions for 
these two runs are shown, respectively, in Figures IV. 8 (A) and (B) . It 
can be seen that the wind direction was fairly constant during the meas- 
urements for both runs. For comparisons, anemometers were also mounted 
on a meteorological tower which was located between two lines of sight 
of the photometers (Figure IV. 6). 

The analog time histories of the a.c. signals are first digitized 
at a rate of 80 samples per second. These digital time histories are 
then averaged every 10 samples to reduce the sampling rate to 8 samples 
per second by the sample reduction method discussed in section IV-3. 

The purpose of the operation is to reduce the quantization errors due to 
the analog/digital conversion. These reduced digital time histories and 
their time derivatives are then used to evaluate the cross-correlations. 
Since the time histories are too long to be stored in their entirety in 
the computer, this evaluation is performed by a specially devised com- 
puter logic for piecewise correlation [IV-14, -15, and -16]. In short, 
a long time history of length T is to break into m shorter pieces of 
equal length AT; T = m AT. Each short piece of data is then processed 
to obtain a "piecewise" cross-correlation. These "piecewise" cross- 
correlations are sequentially averaged by a recursion formula to obtain 
the accumulative cross-correlation which is equal to the cross-correla- 
tion evaluated directly from the entire time history. By so doing, the 
computer needs only to store a time history of length AT. In processing 
the test runs I and II, the piece length AT = 450 seconds, the maximum 
time lag % ax = 50 seconds and the increment in time lag At = 0.25 second. 

Figure IV. 9 shows the cross-correlation of the photometer outputs 
for run I with the lower and upper frequency cut-off, 0.01 Hz and 3.0 Hz, 
respectively. The extremely erratic curve indicates clearly that the 
upper frequency of 3 Hz is too high. An examination of the correspond- 
ing gain shows that a narrow-band noise component exists at 1.25 Hz. So, 
to smooth the cross -correlation, the noise component has to be filtered 
out. Figures IV. 10 through IV. 12 show three different cases of cross- 
correlations Rxy(x) and R^ 2 ) (t) of the outputs and their time derivatives 
of run I, which are obtained by using three different roll-offs in the 
digital filtering. All three have the same terminate frequency frj,=1.0 Hz, 
while the cut-off frequencies, f c , equal 0.8, 0.6 and 0.4 Hz, respectively. 
Also shown in these figures are the 90-percent confidence limits of R X y(r) 

and R^y^(T) as denoted by R m (r) and R^ 2 ^ (x) , respectively [ IV- 15 ] . All 

three R^Ct) curves are much smoother compared to the curve in Figure IV. 9, 
and there is only one peak in R X y(x) around % = 10 seconds, which exceeds 
the confidence limit at least by two times. This means that only this 
peak has the sufficient statistical significance, and that the rest of 
the curve is due to nothing but noises which have not been completely 
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( s') 

eliminated. Therefore, only the portion of (t) which corresponds 

to this significant peak should be selected with sufficient statistical 
confidence. For easy identification, these portions of R.( 2 )(t) are 
shaded. 

Figure IV. 13 presents the statistical errors of of the cross- 
correlations R X y(T) versus the inverse square root of the integrating 
time T = mt. The statistical error is defined as the average of the 

confidence limit R m ( t) integrated over the time lag from -T max to +i max 
[ IV- 14 ] . For statistical stationary processes, the theory shows that 
these curves should approximate a straight line drawn through the origin 
of the coordinates, and that the slope of the straight line gives the 
inverse square root of an equivalent bandwidth B e ff, as if the associated 
spectrum were a white noise spectrum [IV-16]. The figure shows the 
equivalent bandwidth B e ff = 0.182 Hz for run I. It is also noted that 
the statistical error for f c = 0.8 Hz is the smallest, while the error 
for f c = 0.4 Hz is the largest. This is expected, since the error is 
proportional to the inverse square root of the bandwidth used. 

Figures IV. 14 through IV. 16 show the gains G X y(f) and G^y^ (f ) corre- 
sponding to Figures IV. 10 through IV. 12. Table IV-2 summarizes the 
approximate power law for each gain and the corresponding increase in 
exponent from G X y(f) to G^y^ (f ) . For stationary processes, the increase 
in the exponent of the power from G X y(f) should be 2 (see equation (46)). 
Table IV-2 indicates that the case with f c = 0.4 Hz and f T = 1.0 Hz gives 
a power law closest to that for stationary processes. It also implies 
that the components of the photometer outputs between 0.4 Hz and 1.0 Hz 
is less stationary than those between 0.01 Hz and 0.4 Hz. To further 
demonstrate the nonstationarity, G^p) (f ) and f 2 G X y(f) for each case are 
plotted in Figures IV. 17 through IV. 19. Comparing of these figures also 
shows that the best agreement between G^p) (f) and f 2 G X y (f) is for the 
case with f c = 0,4 Hz. 

Finally, we come to check our ultimate goal, i.e., to retrieve 
optically the turbulent wind statistics. Figure IV. 20 shows the normal- 
ized probability density of wind component variations (normal to both 
beams) deduced from the optical cross-beam measurements, as well as 
from the anemometer readings. Clearly, the most probable wind component 
obtained optically is within 5 percent of that obtained by the direct 
measurement. The rms of the wind component variations obtained optically 
is slightly larger than that by the anemometer. The best agreement between 
them again is offered by using f c = 0.4 Hz and f T = 1.0 Hz. It is also 
noted that the probability density obtained by the cross ed-beam technique 
is somewhat skewed to higher wind component compared to that by the ane- 
mometer. The reason for this skewness is still being investigated. 

Figure IV. 21 shows another comparison of the probability densities 
obtained both optically and directly for run II with f^ = 0.01 Hz and 
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f u = 0.3 Hz. Once again, the cross ed-beam measurement yields very 
accurately the most probable wind component. The rms of the wind 
component variations is about 25 percent less than that by the 
anemometer. The optical probability density also shows a similar 
skewness to higher wind components. 


IV- 6. CONCLUSIONS 


The crossed-beam technique employing the natural scattered radia- 
tions is shown to be capable of remote detection of the wind component 
variations normal to the two light beams. A theoretical relationship 
is derived which approximates the probability density of the wind com- 
ponent variations in terms of the measurable cross -correlation of the 
photometer outputs under three fairly realistic assumptions on the 
atmospheric turbulence and wind statistics. Specifically, it is shown 
that the probability density of the wind component is directly propor- 
tional to the temporal cross -correlation of the first time derivatives 
of the light intensity fluctuations multiplied by the time lag intro- 
duced between the two time histories. The moments of the wind component 
variations can be easily evaluated once the probability density is 
obtained by the optical crossed-beam technique. 

Because of the very low s ignal- to-nois e ratio for the time histories 
intrinsic to this type of remote detection, the digital analysis is far 
more preferable to analog analysis for higher statistical accuracy. 
Furthermore, because of the necessity to use the time derivatives of the 
photometer outputs which are approximated by the finite difference scheme, 
the sample reduction method is used to reduce the quantization error by 
about 30 percent. More importantly, a suitable digital filtering should 
be used to filter out the higher frequency noises which are weighted by 
a factor of frequency square due to the necessity of using the time 
derivatives. The digital filtering process is fairly time-consuming 
even for fast computers, and the computation time required is shorter 
for slower roll-off in the frequency response function of the low-pass 
filter. Three different types of roll-offs are studied to see their 
effects on the evaluations of the cross -correlation and the corresponding 
spectral densities. These roll-offs have the same terminate frequency 
of 1.0 Hz, while the cut-off frequencies are 0.8 Hz, 0,6 Hz and 0.4 Hz, 
respectively. 

Actual field tests have been carried out in the atmosphere at 
Boulder, Colorado. Two lines of sight have been crossed about 61 meters 
from the ground. The photometer signals were processed by a special 
piecewise computer program, together with two subroutines for applying 
the sample reduction method and for the digital filtering. The optically 
measured probability densities of wind component variations were found 
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to be in reasonable agreement with those measured directly by the con- 
ventional anemometers mounted on a meteorological tower at the same 
crossing height. The most probable wind components obtained by the 
cross ed-beam technique are within 5 percent of that obtained by the 
anemometer. The root mean squares of the wind component variations 
obtained by the former for two different tests are within 5 and 25 
percent, respectively, of those measured by the latter. For still- 
unknown reasons, the probability density curve obtained by the optical 
means is relatively skewed to higher wind speeds as compared to that 
from the anemometer. It is also found that the best agreement in the 
probability density of wind components was obtained by using the slowest 
roll-off, i.e., f c = 0.4 Hz and f^ = 1.0 Hz. This strongly suggests 
that the best frequency range might be between 0.01 Hz to 0.5 Hz for 
these tests. Theoretical considerations predict that the optimum cut- 
off frequency f c will be determined by two factors: (a) the maximum 

wind component to be expected U max and (b) the smallest atmospheric eddy 
size L^n which will preserve sufficiently its characteristics over the 
transmission between two light beams. Specifically, f c = U max /2jt Lmin- 
This gives f c = 0.318 Hz for U max = 20 m/sec (run I) and L^n = 10 meters 
[IV-8], which is close to the experimentally determined value of 0.4 Hz. 
More experimental results are desirable before we can establish conclu- 
sively an optimum frequency range and the shape of frequency response 
function of the digital filtering with minimum required computation time. 

The analytical results and data-process ing methods developed for 
the atmospheric turbulence measurement can be applied to measure turbu- 
lence in other fluid flow problems as well. The optical technique is 
particularly advantageous over the hot-wire and the like in measuring 
the flow speed and turbulence intensity of high-temperature, chemically 
reacting or supersonic flow fields associated with rocket engine plume, 
boiling phenomena in nuclear reactors, plasma, etc. All of these flow 
fields, at the present time, defy accurate and f ree-from-inter f erence 
measurements by conventional direct-sensing instrumentations. 
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TABLE IV- 1. TEST RUN SPECIFICATION 


SPECIFICATION 

RUN I 

RUN II 

BEAM SEP ACTION 

93 meters 

117 meters 

CROSSING HEIGHT 

51 meters 

61 meters 

ANEMOMETER LEVEL 

62 meters 

61 meters 

BASELINE BETWEEN PHOTOMETERS 

220 meters 

220 meters 

BEAM ELEVATION 

45° 

45° 

BEAM AZIMUTH 

4° to N. 

14° to N. 

MEASURING TIME: 



START 

1:43 p.m. 

11:01 a.m. 

END 

2:36 p.m. 

12:03 a.m. 

SKY CONDITION 

clear 

clear 


TABLE IV- 2. APPROXIMATE POWER LAWS OF GAINS OF 
PHOTOMETER OUTPUTS AND THEIR TIME DERIVATIVES 


DIGITAL FILTER 
SETTING (Hz) 

EXPONENT IN THE POWER LAW 
APPROXIMATION 

INCREASE IN THE EXPONENT 
FROM 

G (f) to G (2) (f) 
xy xy 

f 

c 

f T 

G (f) 
xy 

G (2) (f) 
xy 

0.8 

1.0 

-1.45 

-0.25 

1.20 

0.6 

1.0 

-1.58 

-0.22 

1.36 

0.4 

1.0 

-1.83 

0 

1.83 





Top View 



FIGURE IV. 1. GROUND-BASED CROSSED-BEAM PHOTOMETER ARRANGEMENT 





FIGURE IV. 2. CROSSED-BEAM SIGNAL AND NOISE SOURCES 
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FIGURE IV. 4. ILLUSTRATION OF THE ASSUMPTION 
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FIGURE IV. 5. RESPONSE FUNCTION (PHASE AND GAIN) 
OF THE MARTIN-GRAHAM LOW-PASS DIGITAL FILTER 
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FIGURE IV. 7 


BLOCK DIAGRAM OF THE DETECTION STAGE 
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FIGURE IV. 9. CROSS-CORRELATION AND CONFIDENCE LIMIT FOR 
f, = 0.01 Hz and f = 3.0 Hz (RUN I) 









FIGURE IV. 11. CROSS- CORRELATION AND CONFIDENCE LIMITS 
FOR f_ = 0.6 Hz and f T = 1.0 Hz (RUN I) 




Time delay, T (sec) 


FIGURE IV. 12. CROSS -CORRELATIONS AND CONFIDENCE LIMITS 
FOR f = 0.4 Hz AND L = 1,0 Hz (RUN I) 
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FIGURE IV. 13. STATISTICAL ERROR OF CROSS-CORRELATION (RUN I) 











FIGURE IV. 17. COMPARISON OF G%p (f ) WITH f 2 G x 
FOR f_ = 0.6 Hz AND f T = 1.0 Hz (RUN I) 






0.75 



MM 


CM 

m m 


trj 


H 


H 

Pd 


’ — ✓ 

/—S 

N 

MM 

w 



/ N 

o 

CM >» 

• 

w X 

i — i 

o 

II 

Pm 


o 

H 


m 

s 


o 

Q 

CO 


M 

<3 

ptf 


<3 

N 

PM 

Pd 

§ 


o 

VO 

o 

• 


o 

• 

II 

00 


i — i 

o 

• 

M-l 

> 


M 

Pd 


o 

m 

Pm 

&4 


!=> 


o 


M 


Pm 



104 






FIGURE IV. 20. PROBABILITY DENSITY OF WIND COMPONENT VARIATIONS (RUN I) 
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CHAPTER V 


ON THE TRANSFORMATION BETWEEN TRANSIT DISTANCES 
AND WIND COMPONENT DISTRIBUTIONS 


by 

Fritz R. Krause 

NASA - Marshall Space Flight Center, Alabama 


V- 1 . INTRODUCTION 


The recent conceptual development of the crossed-beam concept indi- 
cates that the correlation of A. C, -coupled photometer outputs can be 
employed to approximate the probability distribution of either transit 
times or transit distances between spatially and temporally separated 
lines of sight. Applications of these measurements then require us to 
relate these transit times or distances to the velocity components assoc- 
iated with the transport of heat, mass, or momentum. The objective of 
this memorandum is to extend this transformation to the use of transit 
distances instead of transit times. 

Present experiments have been restricted to the use of transit 
times. These transit times may be uniquely related to wind components, 
if the direction and length of the transit distance are known. In wind 
tunnel and jet applications, the transit distance is set by separating 
two mutually perpendicular beams along the known flow direction [ V— 1 ] . 
Subsequent transformation from transit times to longitudinal velocity 
variations has successfully matched pitot tube and hot-wire measurements. 
The early field tests of ground winds at tower heights [ V— 2 ] employed an 
outside wind direction reading with a wind vane to obtain the transit 
distance. The wind direction at the top of the tower was measured before 
the test run, and both beams were then rotated normal to this direction. 
The beam separation at the tower height is then a good estimate of the 
transit distance provided that the wind direction fluctuations during the 
run do not exceed ±20 degrees [V-3]. 

Analytical studies of multibeam geometries were initiated to deter- 
mine whether suitable test geometries could be found which do not require 
independent wind direction readings for relating transit times to wind 
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components. W. H. Heybey [V-4] in his investigation of three ground 
telescopes showed that the two horizontal wind components could be 
determined at preselected altitudes, provided that vertical wind com- 
ponents are negligible. A recent extension of his work [V-5] also 
indicates that the third telescope may be replaced with one additional 
photodetector in the focal plane of the second telescope. This simple 
single-beam fan geometry will measure horizontal winds over a restricted 
range of wind azimuths (= ±22°). Vertical wind measurements are beyond 
the capability of both the single-beam fan arrangement and the three- 
telescope arrangement. They would require the use of at least four 
ground-based telescopes. 

The purpose of this paper is to present physical arguments that 
wind components may be uniquely determined with a simple two-telescope 
arrangement if a plurality of photodetectors is used in the focal plane 
of one telescope to monitor transit distances instead of transit times. 
The use of a third or fourth telescope and the restrictions on wind 
azimuth are not necessary, if the measurement objectives are to measure 
one wind component rather than two or three. 


V-2 . THE MEASUREMENT OF TRANSIT DISTANCES 


Consider the test geometry of Figure V.l. One distant radiometer, 
B, is set to look at almost horizontal lines of sight which are all con- 
tained in the same vertical plane. A second radiometer, A, is set 
normal to this plane such that its horizontal line of sight intersects 
the lowest beam of the distant telescope, B. The fan angles between 
the other lines of sight of radiometer B are assumed to be so small that 
deviations from the horizontal may be neglected in the vicinity of the 
intersection point. Each line of sight of radiometer B may thus be 
characterized by its altitude z = h, h 2 , h 3 , ..., above the intersection 
point. Transit distances are measured with the above arrangement by 
differentiating the output signals i^(t) and ig(t,h) of the two radiom- 
eters by delaying signal i^ relative to the signals ig, and by adding 
the lagged product of these time derivatives. The result of this cross- 
correlation calculation is then multiplied with the selected time delay 
and normalized with the resultant peak value as shown in Figure V-2. 
Plotting this normalized multibeam product mean value against the 
abscissa h/x should then directly approximate the probability P(w = h/ t) 
that the vertical speed w of optical disturbances has the value h/x. 

The above claim of approximating vertical wind speed variations is 
based on the same physical model of the common path that was used by 
Heybey. The common signals, which are retrieved by the correlation 
calculation, are identified with optical disturbances which intersect 
beam A at time t - x and B at time t. These disturbances therefore 
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travel along a "common path" that intersects beam A and beam B at the 
altitude h above beam A. Let us assume that the vertical wind component, 
w, does not change along the common path. A disturbance can then be 
"common" to the altitudes z = 0 and z = h only if the vertical transit 
distance component wx is equal to the given altitude h. 


w = h/x. 


( 1 ) 


The time history of horizontal wind component, v, will determine the dis- 
placement of the common path intersection of beam A relative to the beam 
intersection point at z = 0. This displacement distance is given by 

f v dt. 

J 

t-x 


Let us now assume for a moment that the turbulence is homogeneous over 
such distances. The common signals which are generated at this inter- 
section are then statistically equivalent to the disturbances for any 
path intersection point on beam A. In other words, the displacement of 
the path intersection along beam A does not change the product mean value 
calculations. The magnitude of variations of the horizontal wind com- 
ponents v are thus suppressed in the product mean value calculations. 

Any displacement of a common disturbance along the beam cannot be detected 
because of the optical integration along the beam. The above arguments 
may be repeated for the displacement 


t 



t-x 


dx. 


of the common path intersection along beam A at altitude z = h. As a 
result, we thus find that any horizontal motion is cancelled by the 
optical integration along the beams in fields of homogeneous turbulence. 
These same arguments should hold, in first approximation, for inhomo- 
geneous fields of turbulence, where the linear dependence of point-beam 
product mean values on translational disturbances 


v dt 


or 



u d t 


is cancelled by the subsequent integration over all such displacements. 


113 


The above arguments are completely general in that they do not 
depend on the special orientations of the beams B or beam A. The 
arrangement of horizontal beams that is illustrated in Figure V-l was 
chosen to demonstrate that vertical wind components can be obtained 
with only two telescopes. The same arrangement would also closely simu- 
late the desired airplane instrument-package [V-6], However, the partic- 
ular crosscorrelation calculations given in Figure V-2 will indicate 
transit distance components, h, across the fan regardless of orienta- 
tion of the fan normal. For example, horizontal wind components would 
be measured if the beam configuration of Figure V-l is rotated around 
their common baseline until the fan beam normal is in a horizontal 
position. In any case, the crosscorrelation should indicate the proba- 
bility of normal wind components across the fan of beams B 


P(w = h/x) 


const — < 

h 




St 


St 


> 


( 2 ) 


The transformation of these transit distances to wind components is 
already included in equation (2). This transformation no longer 
requires an independent measurement, since the transit time is equal 
to the time delay t = const., which was set by the computer. 


V-3 , ACCEPTABLE RANGE OF TIME DELAYS 


The above discussions would be incomplete without considering the 
restriction on the acceptable range of time delays. A lower limit for 
the time delay, T m i n ) is given by the condition that the associated 
transit distance, wt, should exceed the class interval of its probabil- 
ity distribution. This class interval is given by the length scale of 
the common disturbance, and was reduced to the turbulent microscale, L^, 
by correlating the signal derivatives Bi/St, instead of the signals 
[V- 7] . 


L h 


w 

2jt 


< L 


(t) i B (t; h=0) > 


< 


cjiA(i) 

St St 


> 


( 3 ) 
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The lower limit of transit distances is thus provided by the condition; 


wt . 
mm 


1 

T . 
min 



(t) i B (t; h=0) 
Si A (t) di B (t,0) 
St St 


> 


> 


« 1 . 


(4) 


An upper limit of the time delay, T max , is given by the condition that 
the transit time should be much smaller than the eddy lifetime, L T . A 
first approximation of this eddy lifetime may be derived from the 
envelope to the various time correlation curves, as illustrated in Fig- 
ure V- 3 . 


x § L . 

max x 


Both inequalities together define the range of acceptable time delays; 


2 rt 


i- A (t) i B (t; h=0) > 
Si A (t) Si g (t; 0) 


« T « L 


St 


St 


( 6 ) 


Within this acceptable range of transit times, one should expect that the 
crossed-beam approximation of statistical wind-component distributions, 
equation (2), is valid. One may thus verify this approximation by repeat- 
ing the crosscorrelation calculations for various time delays, x =? t x ; t 2 ; 
t 3 . . . , within the acceptable range. The resulting points should all 
collapse on one curve. The same condition may also be used to obtain a 
sufficient number of points from comparatively few (=5) beams B. , This 
additional advantage is also illustrated in Figure V-2 for two selected 
time delays, x - Tj_ and x = t 2 - The same procedure might be used to 
develop an on-line wind component display system. 


V-4 . ACCEPTABLE RANGE OF FAN ANGLES 


An arrangement of photodiodes in the focal plane of telescope B will 
give a fan of narrow conical fields of view. These fields of view replace 
our hypothetical parallel beams in an actual photodetection system. We 
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therefore have to consider a finite range of fan angles, B, and field of 
view, Af3. The previous idealization of parallel beams may be considered 
as the limiting case, where the fan angles all approach zero. 

The effects of a finite fan angle, p > 0, are illustrated in Fig- 
ure V.4. Consider the horizontal displacement, 


tv u 2 + v 2 , 

across a single photodiode's field of view. The normal displacement, h, 
of a common signal will be subject to a geometrical uncertainty, 


7T T tan 


p s/a 


r, 2 


+ v c 


due to the finite fan angle. This uncertainty should be smaller than the 
tolerance, L^, of the normal displacement distance, h. This gives the 
following geometrical condition for the acceptable range of fan angles, 

0 : 


— T tan p 4 u 2 + v 2 = (7) 

The effects of a finite viewing angle, Ap, will be to generate a finite 
beam diameter, hAp/sin p, which should also be smaller than the toler- 
ance on transit distance. This leads to a geometrical condition for 
the acceptable field of view. 


h AB 

sin p 


V 


( 8 ) 


However, the ratio n = h/L^ gives the number of beams in the fan. Equa- 
tion (8) may thus be rewritten as 


AP 


1 . 

— s m 
n 


0 - 


(9) 
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Equation (7) becomes 


tan 


2L h h = 2 w 

^ T \/u 2 + V 2 11 \/ U 2 + V 2 


( 10 ) 


In case of vertical wind measurements, w / \/ u 2 4- v 2 has values of ^0.1 and 
the largest acceptable fan angle of a five-beam fan would be 


tan p = 



0.1 « 2 °. 


The associated viewing angle would be 
Af3 § | = 0.4°. 


Such a system can obviously be achieved by installing a variable-fiber 
optics block in the already existing detector fan. 


V-5 . CONCLUSIONS 


A plurality of no more than five photodetection elements in the 
focal plane of a single telescope is proposed to monitor wind components 
in preselected directions by monitoring transit distances instead of 
transit times. Physical arguments are presented which lead to the follow- 
ing tentative conclusions. 

(1) Any component of the atmospheric motion along one of the 
beams will not affect the result of crosscorrelation calculations, since 
a displacement of a common disturbance along a beam cannot be detected. 

(2) The crosscorrelation calculations indicate only the com- 
ponents of the transit distances which are normal to both beams. 

(3) The measured transit distances can be directly transformed 
into wind components, since the transit time is known from the time delay, 
which was set in the computer. 

(4) Both horizontal and vertical wind components can be deter- 
mined using no more than two telescopes. 
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(5) The above wind component measurements require us to 
approximate the probability density of wind component variations for 
various time delay parameters. The range of these parameters may be 
determined experimentally with cross ed-beam systems. 

(6) A plurality of photodetectors will always give a finite 
fan angle and a finite viewing angle; whereas, the analysis was based 

on the idealized case of a vanishing fan angle. However, because of the 
finite size of the turbulent microscale, small fan angles and viewing 
angles, which are within the capability of the already existing detector 
fan, are tolerable. 


(a) SID 



119 


FIGURE V.l. CROSSED BEAM TEST ARRANGEMENT FOR TRANSIT DISTANCE MEASUREMENT 




a i A ( t - T ) ai B (t ; h) 



FIGURE V. 2 . PROBABILITY DENSITY OF WIND COMPONENT VARIATIONS 
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FIGURE V. 3 . APPROXIMATION OF EDDY LIFETIME 
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CHAPTER VI 


CROSS- POWER SPECTRAL ANALYSIS FOR ATMOSPHERIC CROSS-BEAM DATA 


by 

J. Briscoe Stephens and John M. Locklin 
NASA - Marshall Space Flight Center, Alabama 

and 

John C. Pooley 

Northrop - Huntsville, Alabama 


VI- 1. INTRODUCTION 


This is the first attempt to physically interpret cross-beam 
atmospheric results using spectral analysis. The cross-beam program 
has used just the time domain for the retrieval of wind velocity infor- 
mation [VI-1, -2], While the spectral formulation applied to this 
analysis is not new [VI-3], the method of combining is new. 

The concept of a piecewise spectrum is introduced, along with the 
statistical error of the spectrum. This paper analyzes the above con- 
cepts and applies them to actual cross-beam data. 

Chapter VII will give the details of how these techniques can be 
applied to retroactive filtering. 


VI- 2 . NOMENCLATURE 

INDEPENDENT VARIABLES 


Symbol 


Definition 

t 

total 

observation time 

T 

total 

integration time 
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Symbol 


Definition 


AT 

length of time for each piece -- the whole 
data record is subdivided into pieces for 
data processing 

i 

piece number 

m 

total number of pieces of data integrated 
over (m = T/aT) 

f 

frequency 

X 

time lag between the two signals 

T m 

maximum time lag 

DEPENDENT VARIABLES 

x(t) 

signal time history for the sensor 

y(t-T) 

delayed time history for the second sensor 

^(x) 

covariance of the ith piece of data at time lag x 

r i (x) 

accumulative covariance average over the first 
ith piece of data 

r(x) 

accumulative covariance average over the entire 
integration time 

At (x) 

statistical error of the covariance 

C(f) 

co-spectra (real part) 

Q(f) 

quad-spectra (imaginary part) 

G(f) 

gain, G 2 = C 2 + Q 2 

<G(f)> 

mean gain over m pieces of data 

G(f) 

accumulative gain over m pieces of data 

<ZG(f) > 

mean statistical error of the gain 

A3(f) 

accumulative statistical error of the gain 
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Symbol 


Definition 


V G,0 the Fourier transpose of the spectra to the 

correlation curve will have the subscripts of 
those spectra components used 

r (t) = r (-t) this is the cross-covariance (subscript xy 
v ‘ indicates positive time lags and yx indicates 

negative time lags 


OPERATORS 

( ) 

A 

~ = S /< > dt 

(i-l)AT 

k 

« >1 

i=l 

“ = <~k=m> 


statistic 


piecewise mean (note: 
real time or time lag, 


accumulative mean where 
data (when average over 

total accumulative mean 


when average done over 
the symbol < > is used 


k is the kth piece of 
pieces ) 


<( ) k > 2 = k iy2 J ( < ‘ ^ i " ( k )2 statistical error 
i=l 

<( )> = F {<( )>] the symbols used in the time domain will be 

transferred into the frequency domain showing 
the origin of the terms 

VI- 3 . FORMULATIONS 
A. Spectra 

Two types of spectrum routines are used in this analysis: the 

accumulative spectra and the mean spectra. The accumulative is obtained 
from a Fourier transformation of the final accumulative correlation curve, 
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whereas the mean spectra is obtained by averaging the Fourier transforma 
tion of each piece correlation curve. The spectrum obtained by these 
formulations is not the common method of spectral analysis. Analog and 
most digital methods of spectral analysis use a fast Fourier transform 
technique. This means that the spectrum is obtained from the raw data 
(before correlation of the data). Since we are Fourier-transforming the 
crosscorrelation curves, we have the spectrum from only the common modu- 
lations of the two signals. Therefore, for large noise, the advantage 
is that most of the noise is cancelled out before the spectrum is 
taken. 


The formulation used to obtain the accumulative cospectra C(f) is 

T 

m 

c(f) = j t r xy( T ) + r yx ( T ) 3 cos dr, (1) 

o 


where x is the delay, r(x) is the accumulative covariance, and f is the 
frequency. The accumulative quad-spectrum Q(f) is 


Q (f ) = / f r xy ( T ) ‘ r yx (T)l sin 2n f x dx. (2) 

o 


It then follows that the modulus of the accumulative gain G(f) is 


GCf) 


J 


r 2 


i , 


( 3 ) 


and the accumulative phase 0(f) is 


0(f) = arctan 


am 

C(f) 


( 4 ) 


These values given by the last four equations are calculated in the com 
puter program given in the appendix. 
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The values for the piecewise gain G-^, co-spectra C^, and quad- 
spectra are calculated in exactly the same manner as their accumula- 
tive counterparts except that the covariance curve used is now the 
piecewise covariance curve. With these terms, the mean gain can be 
defined as 


<G i <£)> 4 I (5: 

The mean gain is a unique term not normally found in the literature. 
There are reasons to believe that the mean gain is free from small 
oscillations in time delay. 

After the error criteria are considered, these two gains will be 
tested on actual data. 


B. Statistical Spectra Error 

One of the key problems in spectral analysis is the lack of a good 
error routine. If one deals with either the co-spectra or the quad- 
spectra, then the standard error routines give unique results. How- 
ever, this is not true for either the gain or the phase, since these 
errors depend on (a) which average (accumulative or mean) is used, and 
(b) the differences in error propagation with changes in integration 
time and summation. 

The statistical error of the accumulative gain AG(f) is defined as 


AG(f) 



( 6 ) 


and it follows then that the statistical error of the mean gain <^G(f)> 
is defined as 


<AG(f)> = 



m 

^[G”(f7 - <GCf7>] 2 . 

i=l 


( 7 ) 


These error calculations were performed using the computer program 
given in the appendix. 
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C. Inverse Transform of Spectra 


The next question is, what happens when the spectrum is transposed 
back into the time domain? This, in effect, does provide a check on 
the work done in the frequency domain. 

Two types of inverse transforms exist: a transpose of the 

co-spectra and quad-spectra, and a transpose of the gain and the phase. 
The component transpose is given by 


f 


max 


r==(T) 

CQ 



f . 
mm 


[C(f ) cos (2jt fx) - Q(f) sin (2jt fx)] dt, 


( 8 ) 


and the gain transpose is given by 


f 

max 



<G(f)> cos [2jt fx + 0(f)] df. 


(9) 


The subscripts indicate the source of the terms used in the transpose. 

The computer program for the above calculations is given in the appendix. 


VI-4. RESULTS 


The formulations discussed in the preceding section will be tested 
on actual atmospheric data from a separated cross-beam run with high 
statistical accuracy to perhaps gain a reasonable empirical understand- 
ing of the spectra. 

The first formulation calculated was the accumulative gain (equa- 
tion (3)) and its associated error (equation (6)). The results are shown 
in the top graph of Figure VI. 1. These results show that the accumulative 
gain is not statistically significant over the range of frequencies from 
0.01 cps to 1.0 cps . We have required only about a 65 percent level of 
probability in the test of significance. Since this run contains a peak 
in the correlation curve with a statistical probability of better than 
99 percent, the spectra from other curves are going to be even less 
significant. To check for errors, the spectrum was transposed back to 
the time domain by equation (8). These results are shown in Figure VI. 2. 
One observes that the transpose of this spectrum gives the initial cor- 
relation curve; therefore, the accumulative spectrum must have the same 
information in it that was contained in the original curve. 
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Before drawing any conclusions we tried another spectral routine. 
Rather than using the accumulative spectrum, we used the mean spectra 
(equation (5)) with its associated error (equation (7)). The results 
are shown in the lower portion of Figure VI. 1. That the mean gain is 
statistically significant suggests that perhaps we are on the right 
track. To verify these results, the mean gain and the mean phase were 
transposed back to the time domain by equation (9). These results, 
shown in Figure VI. 3, appear to be bad, but the mean gain is approxi- 
mately the same as the accumulative gain, according to Figure VI. 1. 
Thus, this result seemed to be only an effect of the variations in 
phase. This theory was tested by transposing both the mean gain and 
accumulative phase, and the accumulative gain and mean phase (see Fig- 
ures VI. 4 and VI. 5). The correlation curve obtained from the mean gain 
and the accumulative phase (Figure VI. 4) gives the best correlation 
curve thus far obtained, and the accumulative gain and the mean phase 
(Figure VI. 5) give about the same results as were obtained for the mean 
gain and the mean phase (Figure VI. 3). 


VI-5. CONCLUSIONS 


Thus far, the spectral results obtained have been analyzed only 
for the most statistically significant atmospheric data we have. Based 
on the empirical results just given, the following tentative conclusions 
are reached. 

The computation of covariance functions from an inverse Fourier 
transform of piecewise-estimated phase and gain spectra may deviate 
from the direct product integration, if the phase spectrum is taken from 
an accumulation of piecewise averages. A reasonable explanation for 
this observation is that the mean gain cancels the effects of nonsta- 
tionarities in transit times that are normally reflected in changes in 
the phase. Thus, evaluating the piecewise correlation curve in the fre- 
quency domain, where shifts in the peak can be neglected, and then 
transposing the final mean result back into time domain may prove to be 
a better technique for dealing with nons tationar ities than just normal- 
izing each piecewise correlation curve. 

Another tentative conclusion is that error criteria can be estab- 
lished only after plotting the error curves for both accumulative and 
mean gain and phase, and comparing the absolute error with the relative 
error. The present error curves should allow the separation of environ- 
mental variations of gain and phase. 
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FIGURE VI. 2. DIRECT AND DOUBLE-TRANSFORMED ACCUMULATIVE CORRELATION CURVES 
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FIGURE VI. 3. DIRECT AND DOUBLE -TRANS FORMED ACCUMULATIVE CORRELATION CURVES 
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FIGURE VI. 4. DIRECT AND DOUBLE -TRANS FORMED ACCUMULATIVE CORRELATION CURVES 
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CHAPTER VII 


CORRELATION CURVES OF RETROACTIVELY FILTERED DATA 


by 

J. Briscoe Stephens and John M. Locklin 
NASA - Marshall Space Flight Center, Alabama 


SUMMARY 


A concept for retroactively filtering a correlation curve is dis- 
cussed and the results given. Normally, the data are filtered before 
being correlated with the net result that, if the limits you have set 
for your high and low pass filters are not suitable, the data must be 
completely reprocessed. Since the data must be either digitized or 
processed through a digital filter and then recorrelated, the time lost 
is considerable. Besides, when the data have been reprocessed, there is 
no assurance that the filters are still correctly set. Further, this 
process tends to make filtering of data in near-real-time wind measure- 
ments of limited value. 

This chapter demonstrates the use of retroactive filtering to pro- 
duce quickly filtered correlation curves from unfiltered data. This is 
done by Fourier transforming the correlation curves into the spectra, 
choosing the frequency range of interest (omitting any instrument noise 
frequencies), and then Fourier-transf orming the spectra back into the 
time domain. 


NOMENCLATURE 


INDEPENDENT VARIABLES 

Symbol Definition 

t total observation time 

T total integration time 
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Symbol 


Definition 


AT 

i 

m 

f 

T 

T 

m 

DEPENDENT VARIABLES 
x(t) 

y(t-r) 

^(t) 

r (t) 

At (t) 

C(f) 

Q(f) 

G(f) 

<G(f)> 

<4G(f)> 


length of time for each piece -- the whole 
data record is subdivided into pieces for 
data processing 

piece number 

total number of pieces of data integrated 
over (m = T/aT) 

frequency 

time lag between the two signals 
maximum time lag 


signal time history for a sensor 

delayed time history for a second sensor 

covariance of the ith piece of data at time lag % 

accumulative covariance average over the first 
ith piece of data 

accumulative covariance average over the entire 
integration time 

is the statistical error of the covariance 

co-spectra (real part) 

quad-spectra (imaginary part) 

gain, G 2 = C 2 + Q 2 

mean gain over m pieces of data 

accumulative gain over m pieces of data 

mean statistical error of the gain 

accumulative statistical error of the gain 
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r G,0 


the Fourier transpose of the spectra to the 
correlation curve will have the subscripts of 
those spectra components used 


yx 


(t) 


= r 

xy 


(-T) 


this is the cross-covariance (subscript xy 
indicates positive time lags and yx indicates 
negative time lags 


OPERATORS 

( ) statistic 



i£T 



(i-l)AT 


dt 


piecewise mean (note; when average done over 
real time or time lag, the symbol < > is used 


« >k> 


< >1 


i=l 


accumulative mean where k is the kth piece of 
data (when average over pieces) 


< > =<( ) k=m > 


total accumulative mean 


k 

<( } k >2 = k^T ^ (( ) t " TT k ) 2 statistical error 
i=l 


<( )> = F {<( )>} the symbols used in the time domain will be 

transferred into the frequency domain showing 
the origin of the terms 

VII. 1 FORMULATION 


The concept of retroactive filtering is straightforward. The corre- 
lation curve is Four ier- transformed into the frequency domain so that the 
unwanted frequencies can be removed. Once these frequencies are removed, 
the spectra are then Fourier-transposed back into the time domain so that 
the correlation curve no longer contains the undesired frequency components 
that the original curve had. 
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The original covariance curve is obtained by 
T 

r xy ( T > = N S " <x> ^f y ( t " r ) ' <y> ] dt - C 1 ) 

o 

This curve is then Four ier-transposed into the frequency domain by con- 
sidering the co-spectra C(f) and quad-spectra Q(f) individually: 


T 

m 

c ( f ) = J' t r x y ( t) + r (t)] cos 2jt fr dr, 
o 


and 


T 

m 

Q (f ) = / [r (t) - r (t)] sin 2jt fr dt. 

x J xy yx 

o 


It follows that the modulus of the gain G(f) is 


( 2 ) 


( 3 ) 


G(f) 



C(f) 2 + Q (f ) 2 , 


and the phase angle is 


( 4 ) 


0(f) = arctan ■ . 

C(f) 


( 5 ) 


The calculations are performed by the computer program given in Appendix I. 
The transpose of the spectrum back into the time domain is given by 


f 

max 


r_ =(t) = 
G 0 



cos 


f . 


mm 


[ 2 * fT + 0(f)] df. 


( 6 ) 
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The frequencies integrated over in equation (6) set the reactive 
filter. For example, if the reactive filter is used to set a band pass, 
only those frequencies within the desired band pass are transposed back 
into the time domain; i.e., f m i n and f max are set at the upper and lower 
limits of the band pass. This means basically that the correlation 
curve arrived at by this practice contains only those frequencies within 
a given band pass. A second form of filtering is sometimes desired; 
notch filtering, which is used to remove a very narrow band pass from 
the data. With the reactive filter, a narrow band pass is removed by 
not integrating over these frequencies. 

The next section will examine how the concept of retroactive filter- 
ing works in practice. This can be performed using the computer program 
given in Appendix II. 


VII-2. RESULTS 


The general procedure of retroactive filtering on a correlation 
curve is summarized from data with a broad band pass (0.01 to 4.0 cps). 

The first case to be considered is the use of retroactive filter- 
ing to set the band pass. The original correlation curve contains' fre- 
quencies up to 4.0 cps. As can be seen in figure VII. la, there is a 
great deal of noise on this curve. The noise frequencies are estimated 
to be between 1.0 and 2.0 cps. Figure VII. lb shows what happens when the 
upper frequency is reduced from 4.0 cps to 2.0 cps. The net result is 
that the correlation curve is almost unchanged, implying that the region 
from 2.0 to 4.0 cps had little effect on the curve. (This curve will be 
used as the basis for the notch filtering.) 

When the upper frequency limit was set to 1.0 cps, the curve in fig- 
ure VII.lc resulted. One can observe how well the retroactive filtering 
has removed the noise on the curve. 

A second example of limiting the band pass of the correlation curve 
by reactive filtering is shown in figure VII. 2. The top correlation curve 
has an upper frequency limit of 4.0 cps, and when the upper limit is 
reduced to 1.0 cps, the bottom correlation curve resulted. It is of 
interest to note how, in this case, reactive filtering actually improved 
the large peak in the correlation curve. 

A second type of filtering commonly used is notch filtering. When 
the data are digitized, certain frequencies are introduced in this pro- 
cess -- if the signal level is low -- that cannot be eliminated by analog 
filtering. Retroactive filtering is a means of eliminating this noise. 
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In the first case (Figure VII. 1) the noise was eliminated by the "brute 
force" method of limiting the upper frequency to just below the digitiza- 
tion noise. The disadvantage with this approach is that the information 
from frequencies greater than 1.0 cps is lost. Thus, ideally we would 
like just to notch-filter the digitization noise. The top curve in Fig- 
ure VII. 3 was notch-filtered in the three digitization noise frequencies; 
i.e., those three frequencies were removed. The bottom curve resulted. 

By increasing the bandwidth of these notch filters from 0.01 cps to 
0.03 cps, the digitization noise can be completely removed. 


VII- 3. CONCLUSIONS 


Retroactive filtering has smoothed the correlation curve. The 
primary advantages to retroactive filtering are that it can be done on 
small computers without reprocessing the data and it saves turn-around 
time. Retroactive filtering is ideal for studies in the regions of the 
spectrum which contribute to transit time measurements. 

The principal disadvantage to reactive filtering is that harmonics 
of the noise will still remain in the results unless they are also 
removed, although in practice the harmonic effect can never be com- 
pletely eliminated. In spite of this disadvantage, it is still believed 
that retroactive filtering is of great value in correlation data analysis 
and can be used in link 3 of the MULTC0R program to make cross-beam 
analysis faster and more effective. 
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FIGURE VII. 3. DIRECT AND DOUBLE -TRANS FORMED ACCUMULATIVE CORRELATION CURVES 




CHAPTER VIII 


THE RELIABILITY COEFFICIENT FOR CROSS- COVARIANCE 


by 

J. Briscoe Stephens and Benjamin C. Hablutzel 
NASA - Marshall Space Flight Center, Alabama 


SUMMARY 


The concept for the reliability coefficient of a covariance curve 
has been developed to improve the previous methods for identifying a 
peak through confidence levels which employ an arbitrarily chosen 
probabil ity . 

The technique described in this Chapter differs from the technique 
traditionally used for peak identification in that, rather than setting 
an arbitrary probability, the probability is now determined as a function 
of the statistical error and the covariance coefficient. For lack of a 
better term, this function has been defined as the reliability coefficient. 

The physical significance and the reliability coefficient are dis- 
cussed in terms of a stationary model. Results from the application of 
the reliability concept to actual atmospheric cross-beam data are 
discussed. 


NOMENCLATURE 


INDEPENDENT VARIABLES 


Symbol 


Definition 


t 


total observation time 


T 


total integration time 


AT 


length of time for each piece -- the whole data 
record is subdivided into pieces for data pro- 
cessing . 
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Symbol 

i 

m 


f 

T 


Definition 
piece number 

total number of pieces of data integrated over 
(m = T/AT) 

frequency 

time lag between the two signals 
maximum time lag 


DEPENDENT VARIABLES 
x(t) 

y(t-r) 

r i (t) 

r (t) 

Ar (t) 

P(tp,m,T) 

OPERATORS 
( ) 

iAT 

/ < > dt 
(i-l)AT 


signal time history for the sensor 

delay time history for the second sensor 

covariance coefficient of the ith piece of data 
at time lag t 

accumulative covariance average over the entire 
integration time 

statistical error of the covariance 
reliability for a time delay 


s tatis tic 


piecewise mean (Note; when average done over 
real time or lag, the symbol < > is used.) 


« >k> " k 


( ). 


accumulative mean where k is the kth piece of 
data (when averaged over pieces) 


()=<() 


k=m 


> 


total accumulative mean 
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Symbol 


Definition 


k 

<( )k> 2 = j~y ) (( )^~( )^) 2 statistical error. 
i=l 


VIII- 1. INTRODUCTION 


The concept of a reliability coefficient for the evaluation of 
cross-covariance curves obtained from atmospheric cross-beam measure- 
ments is developed in this chapter. The reliability concept is a new 
way of viewing the traditional peak identification procedure, which 
uses an arbitrary probability level based on test results. The pro- 
cedure described here requires that the maximum positive peak in the 
covariance curve (the mean covariance) be greater than the statistical 
error. When this stipulation is fulfilled, the time lag corresponding 
to this peak is defined as the transit time. This is the lapsed time 
used by the turbulence eddy in traversing the distance between the beams 

When an arbitrary probability level is used, validity of the corre- 
lation curves is obtained as a binary response and not as a relative 
measurement. Comparison of the statistical error at different time lags 
does not determine which time lag is more statistically significant. 

This fact means that the current approach does not afford a systematic 
means to compare different parts of the same correlation curve. Further 
more, the variability of the experimental parameters from test to test 
prohibits the use of marginal correlation results since the use of an 
arbitrary probability level requires an extra safety factor to account 
for all test conditions. 

The reliability coefficient being proposed here affords a relative 
evaluation of the reliability of the correlation curve based on both the 
sample population and the statistical error. This relative evaluation 
results in the possible elimination of the limitations of the arbitrary 
probability levels when the reliability coefficient is used as the key 
criterion for evaluating the consistency of the correlation results. 

The same information can be obtained directly from the reliability 
probability used with arbitrary probability levels. 

It is hoped that this approach can account for piecewise variations 
in both the transit time and the amplitude of the peak. It should be 
recognized that this is merely a proposed approach and not a solution 
to these problems. 
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This chapter assumes that the reader is familiar with both the corre- 
lation technique (Jayroe and Su [VIII-1]) and the atmospheric cross-beam 
program (Stephens, Sandborn, and Montgomery [VIII-2]). 

VIII-2. THEORETICAL DEVELOPMENT 


The statistical error of the covariance can be interpreted physi- 
cally, as follows, rather than mathematically. 

The accumulative covariance is defined as 


r(T) 



x] [y (t-x) 


y] dt, 


( 1 ) 


where x and y are the instantaneous fluctuations of the signal, x and y 
are the accumulative mean values, T is the total integration time, t is 
the real time measurement from the start of the record, and t is the 
time delay. 

The statistical error At(t) for the covariance is 

W 2 

[r.(x) - r(T) m ] £ j , (2) 

i=l 


m 


£r(x) = 


(m - 1) 


where m is the number of pieces of data used, r^(x) is the ith covariance 
occurring at the time delay x, and r (x) is the resulting accumulative 
covariance for the mth piece of data. 

When using the statistical error for providing a confidence limit, 
it is required that 


r(x) > Ar(x) • t p=g0% 


( 3 ) 


for the covariance coefficient to be significant. This 
807 o, t _n n7 , implies a probability level of 80 percent. 

p OU/o 


Student's "t" of 
There are three 
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categories into one of which a covariance coefficient can fall by the 
scheme set forth by this last restriction. The covariance coefficient 
at a time lag can be less than, greater than, or equal to the statis- 
tical error. 

When the statistical error is greater than the covariance coef- 
ficient, i. e. , 


r (t ) < £r(r ) 

cl d 


*>=80’ 


(4) 


then we say that this is not an acceptable measurement. However, a 
measurement is statistically significant, i.e., 


r(T b ) >Ar(r b ) • t p=8Q , (5) 


if the average ( ) exceeds its own standard error A( ). 

The third category is the one where the statistical error is equal 
to the covariance coefficient: 


r(T c ) = At(t c ) • t p=g0 . (6) 

This relationship implies that there exists an 80 percent probability 
that the sample and the expected value of the covariance have the same 
s ign. 


Thus, when the covariance coefficient of a positive covariance 
sample equals the statistical error, then we have a relation for the 
probability of having a positive correlation coefficient. This rela- 
tion is 


r (t) = At(t) • t 


(7) 


where the probability level, p, is the unknown. However, this unknown 
is related to the known confidence factor 


t = - r -^~ 

P Ar (x) 


( 8 ) 


by the Student's "t" distribution. 
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The probability level required to satisfy the Student"s "t" rela- 
tionship given in equation (8) is 


x=t 


P (m; t) = 


T(m/2) 


n/ Tt(m-l) 


1 1 2 ; x=0 


[1 + 


m 

2 


(m-1) ' 


dx. 


(9) 


The result for P(tp,m,x) will be defined as the "reliability coefficient." 
Although the reliability probability was developed for the cross-covar- 
iance, it is equally applicable to piecewise derivations of any other 
statistics. The Student's "t" distribution was used because the samples 
under consideration are normally small (m < 30). For larger piece 
numbers, the normal distribution might be used, i.e., the reliability 
coefficient becomes independent of the number of pieces used. 


VIII-3 . APPLICATION OF THE RELIABILITY 
TO THE STATIONARY COVARIANCE CURVE 


The variation of the reliability coefficient with time lag can be 
used to afford a more meaningful peak identification technique because: 


(1) 

The reliability curve provides a systematic procedure 
to establish the relative reliability of different 
portions of the same covariance curve. 

(2) 

The reliability curve can be used for estimating the 
noise level tolerance. 

(3) 

The reliability coefficient accurately accounts for 
the uncertainty associated with small sample size by 
employing the Student's "t" distribution. 

Before discussing the properties of the reliability coefficient, 
we must first review additional characteristics of both the covarianc 
and standard error. 


The limitation in computer storage prohibits the integration of the 
covariance (equation (1)) over the entire record at one time. Therefore, 
the record is subdivided into m pieces at length AT such that 


T = m • AT. 


( 10 ) 
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For each piece, a piecewise covariance coefficient r(x)^, where i is 
the piece number, is obtained. Then the accumulative covariance coef- 
ficient r(x) m for the entire record is 


r ( t) = <r ( t)> 

m m 


m 



i=l 


( 11 ) 


Thus, the statistical error Ar(x) defined in equat ion (1) is the effec- 
tive measure of how the piecewise covariance r ( t) £ varies about the 

accumulative covariance r(i) m . The primary physical requirement of the 
correlation technique is that the phenomenon being measured is recurring. 
Thus, the consistency of certain piecewise covariance coefficients is 
the necessary property (implied by the requirements of the correlation 
technique) that enables us to retrieve temporal information. The 
reliability is nothing more than an estimate of this consistency. 

To expand on the significance and application of this requirement 
of consistency, let us consider the following argument. This argument 
will be restricted to stationary data obtained from the convection of 
a frozen pattern at a constant speed such that the eddies traverse the 
transit distance between the two beams in a transit time, t s . The 
piecewise covariance coefficient at the time delay t s should be 


0 < r (T s ) i < 1, 


( 12 ) 


and will be referred to as "the co variance due to the signal." The 
remaining covariance coefficients, t'(T s )£ (the prime indicates all 
time delays other than when t = t s ), can have values between 


“ 1 < r' ( t s )£ < 1, (13) 


and will be referred to as "the covariance due to the noise." (For ease 
of discussion, the signal component is defined as a spike in the covar- 
iance curve. This does not restrict the effect of this discussion.) 
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For random stationary data, the signal component of the accumu- 
lative covariance curve has the following property: 


0 < lim r(T s ) = constant g 1, 

m — >«> 


(14) 


while the noise component has the property that 


lim r ' (x g ) ~ 0. (15) 

m — »<x> 


This property has the resultant implication that the peak in the covar- 
iance curve due to the signal remains at approximately the same time 
lag from piece to piece, but a peak in the piecewise covariance due to 
noise will be moving from piece to piece. Previous peak identification 
efforts thus assumed that the peak corresponding to the convection speed 
in the accumulative covariance curve can be differential from a peak of 
equal magnitude due to noise by using the statistical error of the 
covariance, since the statistical error of the signal peak is less than 
the statistical error due to the noise peak. A more sophisticated 
approach is required in practice because of finite integration times 
and because of the varying magnitudes of the peaks in the covariance 
curve. 

The reliability concept affords an approach which should circumvent 
these problems. To understand the behavior of the reliability coeffi- 
cient, let us consider the following (see eqs . (7) and (8)): A large 

reliability coefficient at a time delay means that the covariance coef- 
ficient has the same sign from piece to piece; whereas, a small reliabil- 
ity coefficient means that the covariance coefficient at that time delay 
is oscillating from piece to piece. These oscillations are the charac- 
teristics of noise and it is thus apparent that the signal should have 
a higher reliability than the noise, once a sufficient sample of the 
population has been obtained. Thus, the reliability curve affords a 
method that, along with the covariance curve, can permit the signal to 
be distinguished from the noise in the curve before the noise has been 
completely suppressed. 
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The above discussion can be extended to set a noise tolerance by 
using a knowledge of the beam geometry and the wind direction. With 
this knowledge, one can determine whether the transit time should be 
a positive or negative time lag. The signal half-plane will be defined 
as the half-plane of the covariance curve in which the transit time 
must exist. The other half-plane is the noise half-plane, since any 
peak existing in this half-plane will be assumed to be due to noise. 

The maximum reliability coefficient \<rhich exists in the noise half- 
plane is then a measure of the minimum acceptable confidence level, 
and will thus give the noise tolerance. 

This approach obviously has some shortcomings, the major one being 
a periodic signal. For a periodic signal, the noise half-plane will 
have a reliability equal to the sum of the noise and the signal. Thus, 
in this case, a more complex formula will have to be used to obtain the 
noise level. An example of such a formula is to subtract the known 
signal reliability coefficient from the maximum reliability coefficient 
existing in the noise half-plane. The result is the noise reliability 
level. These problems should not exist with a broadband process, but 
it is always necessary to know the kind of process that is being meas- 
ured. The above concepts are thus incomplete and should be taken as 
a basis for a new systematic approach to peak identification, not as 
a finished identification criterion. 

To summarize the application of the reliability coefficient for 
setting the noise tolerance, let us consider the following two applica- 
tions. First, for zero beam separation (Figure VIII. 1), the procedure 
for the interpretation is as follows. The signal is what occurs around 
zero time delay. All other peaks then must be noise peaks. If the 
noise reliability coefficient is less (see peak (a)) than the signal 
reliability coefficient, then the covariance curve can be meaningfully 
evaluated; but if the noise reliability coefficient is greater (see 
peak (b)) than the signal reliability coefficient, then it does not 
seem to be meaningful to accept these results. (The case of a periodic 
signal has been omitted.) 

The second application is to the measurement of a wind speed (Fig- 
ure VIII. 2) . In this case, the signal range of time lags is set by the 
range of possible wind speeds during the period of observation. Once 
this is established, the same argument that was used in the case of a 
zero beam separation applies here. If the reliability coefficient of 
the signal is not greater than the reliability coefficient of the noise, 
then the results are questionable. 

The reliability affords one final feature that can be useful for 
setting the integration times. When the reliability of the signal peak 
does not increase (or if it decreases) with additional integration time, 
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and if the probability of noise does not decrease, then additional 
integration time will afford no improvement in the covariance curve. 
The advantage of this over the accumulative error currently used is 
that it provides closer control. That is, the accumulative error 
monitors only the overall variation, while the reliability can be used 
to monitor specific peaks and to determine how the integration time 
affects these peaks. 


VIII-4. APPLICATION OF THE PROBABILITY 
TO ACTUAL ATMOSPHERIC DATA 


In section VIII-3 the discussion was limited to stationary data, 
because at best this approach leaves a certain amount to conjecture. 

When dealing with nonstationary data, the amount of conjecture increases 
therefore, the power of these concepts will be empirically demonstrated. 

The first case shown in Figure VIII. 3 is a model example of the 
use of this probability concept. The peak in the covariance curve at 
-10 seconds compares favorably with the wind speed. The signal peak 
of the reliability coefficient is about 99 percent, while the maximum 
noise peak is only 73 percent. The fact that the reliability coeffi- 
cient of the signal is greater than the reliability coefficient of the 
noise indicates that the signal peak is significant, a fact already 
obvious from the statistical error. This is the only wind measurement 
evaluated that had a 99 percent reliability coefficient for the signal. 
We have found by independent experiments of other statistical techni- 
ques that this is our most well-behaved run. 

The second case of a wind measurement, shown in Figure VIII. 4, 
demonstrates the use of the "noise" half-plane. Winds compare only at 
negative lags. The peak at -15 seconds corresponds to the wind speed, 
but the peaks at -25 seconds is also statistically significant, accord- 
ing to the statistical error. Our concept holds that the noise peak 
at 50 seconds has a reliability coefficient of 71.6 percent. This 
eliminates the peak in the signal plane at -25 seconds, since it has 
a reliability coefficient of only 70.2 percent. We are left with the 
signal peak of a 81.6 percent probability, which is what we desired. 

The third example is a zero beam separation (Figure VIII. 5). 
According to the statistical error, we should reject the results shown 
in Figure VIII. 5, even though they are the desired ones, because they 
are not statistically significant according to the statistical error. 
When the reliability coefficient is applied, it is found that it peaks 
at zero time delay (out to ± 75 seconds). Thus, the reliability coef- 
ficient identifies the peak at zero as meaningful. However, this 
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identification also indicates that the probability for the identified 
peak to be positive is only 21 percent. Our confidence in the "accepted" 
result is thus still low. 


VIII-5 . CONCLUSIONS 


At this time, only a preliminary conclusion is possible since a 
great deal more work is needed with the reliability coefficient before 
a conclusion can be drawn. At the time this concept was tested, the 
concepts discussed in section VIII-3 seemed to be valid for all of the 
cross-beam data available at that time. Furthermore, this concept 
seems to increase our ability to evaluate the covariance curve, and 
seems to work for both stationary and nonstationary data. Tests sug- 
gest that it also may be useful in evaluating derivatives and spectra. 
More work is planned with the latter. 
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Reliability Coefficient 



FIGURE VIII. 1. MODEL RELIABILITY CURVE 
FOR ZERO BEAM SEPARATION 


160 


Reliability Coefficient 



FIGURE VIII. 2. MODEL RELIABILITY FOR A BEAM SEPARATION 
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♦ Statistical Error, 
Covariance Reliability Coefficient 
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FIGURE VIII. 4. ACTUAL BEAM SEPARATION 




Covariance Statistical Error, 

Reliability Coefficient 

.10 --.50 



FIGURE VIII. 5. ACTUAL ZERO BEAM SEPARATION 
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APPENDIX I 


SPECTRA PROGRAM 


This program is designed to calculate the cross-beam power spectral 
density function from the covariance function by using the Fourier 
transform which converts from a time domain to a frequency domain. 

This program was modified from the current spectral program used 
in CORFUN in the following ways: 

(1) It has the ability to shift the major peak to zero. 

(2) It is capable of using Parzen's weighting factor to pre- 
whiten the covariance function. 

(3) It has the option to add zero correlations to increased 
time lag range which provides for a finer frequency 
resolution (1(2* T max ) = Af ) . 

(4) It can provide calculations from covariance coefficients 
taken directly from the output of the CORFUN program, 
either from piecewise data or from the accumulative 
averages of these pieces. 

Table 1 defines the input parameters, maximum run time, and 
available computers at MSFC. Table 2 shows the loading procedure 
for accumulative data, and Table 3 shows the loading procedure for 
piecewise data. Table 4 shows all of the output parameters for 
this program. 


167 



168 




TABLE 1 


COMPUTER: 


SDS 930 - LIMITED STORAGE 

CDC 3200 



CAPABILITY 


RUN TIME / PIECE: 


10 MINUTES 

6 MINUTES 

DEFINITIONS OF 




INPUT DATA: 




INPUT CARD 

© 

RUN 

-RUN NUMBER 


PHOTO ^ 

-1 PHOTOMETER 



PHOTO 2 

-/IDENTIFICATION 

INPUT CARD 

© 

DELT 

-TIME LAG BETWEEN POINTS 


NLAGS 

- NUMBER OF POSITIVE TIME LAG POINTS 



LSCALE 

--IF EQUAL TO 1, MULTIPLIES 
EACH CORRELATION POINT BY 
RSCALE; IF EQUAL TO 0, THIS 
PART IS OMITTED. 



RSCALE 

- A SCALE FACTOR WHICH 
SCALES THE VALUES OF 
THE CORRELATION POINTS; 
USUALLY SET EQUAL TO 1. 



LPARZN 

- IF EQUAL TO 1, CALLS 


PARZEN SUBROUTINE, IF 
EQUAL TO ZERO, IT IS 
OMITTED. 
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TABLE 1 -con't 


INPUT CARD (2) HAGS ENABLES PROGRAMMER TO 

w ADD ADDITIONAL TIME 

DELAYS TO PROVIDE FOR 
A FINER FREQUENCY 
RESOLUTION. 

TAUSHF VALUE WHICH THE PEAK 

IS TO BE SHIFTED. USED 
ONLY WITH PARZENS. 

NORMALLY NOT USED IN 
OUR WORK. USUALLY SET 
EQUAL TO 0. 

LYPE IF EQUAL TO 1, READS 

ACCUMULATIVE DATA; IF 
EQUAL TO 0, READS 
PEICEWISE DATA. 

INPUT DATA © TLAGS REAL ARRAY CONTAINING 

TIME LAGS. 

RXY ARRAY OF POSITIVE 

CORRELATION COEFFIEIENTS. 

INPUT CARD © NTEST IF EQUAL TO - i IT READS 

ANOTHER SET OF ACCUMULATIVE 
DATA AND CALCULATE SPECTRA. 
IF EQUAL TO 0 TAKES 
ANOTHER SET OF PIECEWISE 
DATA. IF EQUAL TO +1, 
TERMINATES. 
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TABLE 4 


PRINTED: 


(TLAGS) 

(RXY) 

(RYX) 



THIS 

IS FOR 

CHECKING 

xxx.xxxxx 

xxx. xxxxx 

l 

xxx .X xxxx 
1 


> 

THE 

ORIGINAL 

INPUT 

1 

1 

1 

1 

1 



DATA 



xxx.xxxxx 

xxx.xxxxx 

xxx.xxxxx 

7 






(PARZEN) 

(RXY + RYX) 

(RXY - RYX) 

\ 



xx.xxxxxx 

1 

xx.xxxxxx 

1 

xx.xxxxxx 

1 


1 THIS OUTPUT 

PERMITS 

1 

1 

1 

1 

1 

1 


> CHECKING OF 

CALCULATION 

xx.xxxxxx 

xx.xxxxxx 

xx.xxxxxx 

) 




FREQUENCY 

COSPECTRA 

QUADSPECTRA 

GAIN 

CROSS PSD 

PHASE 

LOG FREQ. 

LOG GAIN 

LOG CROSS PSD 

(CPS) 





(DEG) 




xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

1 

1 

xxx.xxxxxx 

1 

1 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

xxx.xxxxxx 

xxx.xxxxxx 

1 

xxx.xxxxxx 

1 

xxx.xxxxxx 

I 

xxx.xxxxxx 





END OF PRINTED 

OUTPUT 












PUNCHED 

ON EACH CARD: 







FREQUENCY 

COSPECTRA 

QUADSPECTRA 

GAIN 

CROSS PSD 

PHASE 

LOG FREQ. 

LOG GAIN 

LOG CROSS PSD 

(CPS) 





(RADIANS) 




XX. xxxxx 

xx.xxxxxx 

XX. xxxxx 

XX. xxxxx 

XX. xxxxx 

XX. xxxxx 

XX. xxxxx 

XX. xxxxx 

XX. xxxxx 


OUTPUTS FOR SPECTRA PROGRAM 



SPECTRA PROGRAM 

: SPECT COMPUTES THE cross spectral density function 

equivalence iryx.rmh 

FQLIIVAIENCE (RXY,RP| 1 

DIMENSION RXYM06I , R Y X C 4 0 6 I > TL AGS I 4 0 6 I , PPL t 4 0 6 J , RM I [406 J,C[406I, 
60 r 4 nc, i 

RF AD 1 nj, RUN. PH0T01 , PHOTO? 
mi format rxEio.ni 

READ 3 • DFI-T, Nl A GS , L SC ALE , Li. AGS , L PAR7N, RSCALE , TAIJSHF. LYPF 
3 FORMAT fFI 0 . 0 ,4 TX.Flf) , 0 ,F1 f) . 0 , 13 ] 

C AS THF FOURIER TRANSFORM OF THE CORRELATION 

C COEFFICIENTS. WHERE 

. C 
C 

C ** r NPUT * * 

C N LACS = THE NUMRER OF CORRELATION LAG POINTS, INTFGFR . 

C DEL.T - THE SAMPIING INTERVAL! DELTA TIME!, REAL. 

C Tl AGS = A RE A I ARRAY CONTAINING THE TIMF LAGS. 

IF 1 LYRE 1 38, ,59.38 

38 CONTI NUF 
4 ft n 7 CONTINUE 

00 1 1-1 , Mi AGS 

READ 7 , TLAGSr I 1 , RXYt I J 
? FORMAT MX, ?E1 0.01 

TL AGS r I 1 =ABSF f TL AGS I III 
1 CONTINUE 

no 13 1 = 1 , NL. AGS 
READ 14, RYX ( I 1 
1 4 FORMAT t 1 1 X , F 1 0 . 0 ] 

13 CONTINUE 
GO TO 6 R 

39 CONTINUE 

DO 3/ 1=1, Nl. AGS 

RFAD 46, RXYf I 1 .RYXT I I, T LAGS! I I 
46 FORMAT M X, F10 . 4 , 2E10 . 6 I 
3; CONTINUE 

c RXY = RE Al ARRAY IF POSITIVE CORRELATION COEFFICIENTS. 

C RYX = RE A I ARRAY IF NEGATIVE CORRELATION COEFFICIENTS. 

C 

c **nuTP!IT** 

C R P I ~ a REAL ARRAY CONTAINING THE EVEN PARTt+I OF THE 

C CROSS CORRELATION FUNCTION. 

C RM! - A REAL ARRAY CONTAINING THE ODD PART [-1 OF THE FUNCTION 

C C = RE A I ARRAY CONTAINING CO-SPECTRA, 

C Q = REAL ARRAY CONTAINING QU AD - SPECTRA . WHERE- - 

C G t K I = C r K 1 - J * Q f K 1 

C 

65 CONTI NUF 

1 F ( l I AGS-NI, AGS 1 36 ,36,35 
35 L=Nt AGS+1 

DO 3? I = L , LL AGS 
RX v [ t 1 = 0 . n 
RYX [ I MO . n 
IF II SCALE] 61 ,6?, 61 
6 ? Tl, A G S’ f I 1 s I *DEI T 
GO TO 3 7 
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61 TLAGS [II = ! 

32 CONTINUE 
NLAGS=LLAGS 
36 CONTlNtlF 

PRINT 102, RUN,PH0T01,PH0T02 

102 FORMAT (1 X , 4HRUN=F10 .6, L2HPH0T0MFTFRi=F10.6, 12HPHOTOMETER2=F10 .6) 
PRINT 22, NLAGS, DELT 

2? FORMAT II X,6HNIAGS=T5,5X,5HDFLT=F10.5J 
-PRINT 100 

100 FORMATUX,! OHTIMR L A GS , 4 X , 3HRX Y , 7X , 3HR YX 1 
PI =3 . 14159 
OT = OFL.T 
M = NLAGS 

R IS THF USUAL VARIABLE OF SUMMATION WHICH NORMALLY RUNS R=0,1,... 

M--WHERE M IS THE NUMBER OF LAGS, FOR CORFUN, TAUMIN IS NOT 

NFCFSSARILY ZFRO so THAT THE RATIO R/M RUNS FROM RMIN TO 1.0 
RMIN=TL AGSfl J/TLAGSIMJ 

RMAX = 1.0 
UV = M-1 

RDEL = [RMAX-RMIN1 /UV 

COMPUTE R + AND R- BOTH WFIGHFD BY LAG FUNCTION 

R r RMIN 

IF [TAUSHF] 400,400,401 

401 CONTI NIJF 
MM=TAUSHF/DELT+1. 0 
NLAGS=NL AGS-MM+1 
MsNLAGS 

no 407 1 =1 .NLAGS 
N = I — 1 
MK = N + MM 

RYX IMK1 =RYXf I 1 
407 CONTINIIF 

DO 409 1 =2, M 
J= I -1 
J = .I-MM 

IF r .1 + 1-1 402,40 3,404 

402 J= ABSF I J 1 

RYX I J I =RXY f I 1 
GO TO 409 

403 RXYfl 1 =RX Y f MM 1 
RYX M. I = RX Y I MM 1 
GO TO 409 

404 J=J+2 

RXYI J]=RXY( I I 
409 CONTINUE 
KK = 0 

NN=Ni AGS-MM 
KK=KK+NL AGS 

DO 411. 1 = 1 , MM . ... 

NN=NN+1 

RXY [NN1 =RXY IKK ] 

KK=KK+1 

411 CONTINUE 
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400 CONTINUE 

DO 11 1=1 , NLAGS 

. IF ELSCALEJ 51., 52, 51 

51 RXY t I 1 = R X Y [ T I * R S C A I E 

RYXU I =R YX. I I 1*RSCAIE ,/ 

52 CONTINUE 

PRINT 30, TL AGS Ml. RXY m.RYXm 
30 FORMAT MX.3F10.5I 

11 CONTINUE 

DO 5 I = 1 , M 

C COMPIJTF PAR7EN LAG WEIGHING FUNCTION 

IF1LPAR7N1 40,40,41 

40 D=1 . 0 

GO TO 42 

. 41 CALI PARZEN (R.DI .... 

42 CONT I NUE 
D = 0/2 . 

1JV = RXY til 
T = UV + RYX f I 1 
S = !JV-RYXfIl 
RPL HI = D * T 

rm u 1 1 = n*s 

PRINT A, D , T , S 

6 FORMAT f 3 X , 3 F 1 5 . 8 I 
R = R+RDEI 

5 CONTINUE 

C COMPUTF VARIOUS CONSTANTS FOR TRANSFORM 

TWODT = 2 . *DT 
FORDT = DT *4 . 

. TWOP I = 2 . *P 1 
R = RM 1 N 
N = M-1 
DO 20 1=1 , M 
UV = R *P I 
CARG = UV 
CSUM = 0. 

OSIJM = 0, 

no 1(1 J = 2 , N 
ARG,= CARG 

JFIARG-TWOPJ] 7,7,12 

12 CONTINUE 

IM = ARG/TWOP1 
FIM = IM 

ARG = CARG-F1M*TW0P? 

7 CONTINUE 

CSUM = CSUM + RPI. t J 1 *COSF I ARG ] 

QSUM = OSUM + RMI IJ1*SINFtARG] 

CARG = CARG + tJV 
10 COMTINUF 

0 f II = F ORDT *OSUM 

Cl ! ]=TWODT*IRPL,tll + [ ( - 1 )**[ 1+1 )*RPLfMl ] 1 + E 0 R D T * C S U M 
R = R + RDFI. 

20 CONTINUE 
U V = N I AGS 

IF 'I SC A I F ] 8 , V , 8 



8 TLAGS f M 1 =TLAGS [ H ] *DELT 
TLAGSM 1 = T L A G S f 1 I * DEL T 

9 CONTI NUF 

TML A GtJ = T I AGSIM1 
TMLAGLsTI AGSI1 I 
FMAXs ,5*UV/ tTML AGU-TMI AGI ] 

F M I N = 0 . 

FDEL = f FMAX-FMI N 1 / IJ V 
PRINT 103 

1 0 3 FORMAT f5X,9HFRFQUFNCY,3X,9HC0SPFCTRA,lX,llH0UADSPECTRA>6X. 4HGAIN. 
1 AX . 1 OHCROSS PSD,3X,5HPHASE,4X,9HL0G FREO , 5X . 9HL0G GAIN.3X, 

1 1 SHI. OR CROSS PSm 
DO 30 j 0 I =1».NL AGS 
GAlM = SORTFtC.f I i*cm + OllUQim 
IF ICtI ] 13015.3012,3015 
301 5 CONTINUE 

ph asf = o m/cm 
phasf = at ANF 1 PHASE 1 
GO TO 3013 
3 012 PHASE -1 . 5708 
3013 CONTINUE 

IF [CM ] ! 3333.3334,3334 
3333, PHASF-PHASE + 3 . 14159 

3334 IF r PHASE 1 3 335,3336,3336 

3335 PHASE=PHASE+6 . 28318 
M36 CONTINUE 

PHA=57 ,2958?*PHASE 
CPSD=GA I N*FM IN 
I F r FM I Nl 1 1 001,1001,1 002 
100? CONTINUF 

CLOG-Al OG f CPSDI *0 . 4343 
FI OG = A10G I FM IN 1 +0 . 4343 
1001 CONTINUE 

Gl OG = A I. OG I GA I Nl *0 . 4343 

PRINT 2J . .FMIN.Cfll ,0(1) ,GAIN,CPSD,PHA,FI 0G,6L0G,CL0G 
21 FORMAT r 1 X , 9F1 2 . 6 1 
cr I !=Cf I 1/10 
on l = o r t i/io 

FLOG=FLOG/10 
Gl OG = GI OG/1 0 
CLOG=CLOG/1 0 
23 FORMAT MX.9F8.5I 

PUNCH ?3,FMIN,Cf I I -Of I I , GAIN, CP SD, PHASE, FLOG, GLOG, CLOG 
FM ! N-FM I N + FnEI. 

3010 CONTINUF 

re a n ga, mtest 

56 FORMAT M3) 

IF I M T F S T 1 4 0 0 7,39,4006 
4006 CONTINUF 
STOP 
END 

SURROUT I NF P A R 7 F N f R , D 1 
RC = R 
Q = 1 . -RC 
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IF [RC-.511 0,10,20 
in D = 1 - - f 6 . *RC*RC*Q 1 
RETURN 

?n o = ?. *q*o*o 

RFTURN 

END 
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APPENDIX II 


PROGRAM FOR FOURIER TRANSPOSE 


This program is designed to calculate the covariance function from 
the spectra. This is done by using a Fourier transpose of the gain and 
phase to convert from a frequency domain to a time domain. 

One property of this program is that all or any group of frequencies 
(from 0.01 - 4.0 cps) can be either omitted or used in calculating the 
covariance curve. The inputs must include the beginning frequency and 
the ending frequency along with any certain frequencies to be omitted 
from the calculations. 

Table 5 shows available computers, maximum run time, and a defini- 
tion of input terms. Also shown are the output parameters for both 
printed and punched outputs. Table 6 shows a loading procedure for 
data along with the input parameters. 
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TABLE 5 


COMPUTER : 



SDS 

930 

CDC 3200 

RUN TIME: 



DEPENDENT ON THE FREQUENCY RANGE 
DESIRED-- FROM 0.01-4.0 cps MAX. 
TIME IS 25 MINUTES. TIME IS 
DECREASED BY USING CDC 3200. 

DEFINITIONS 

INPUT 

OF 

DATA: 





INPUT 

CARD 

© 

NLAGS- - 

— 

NUMBER OF TIME LAG POINTS 

WHICH YOU INTEND TO CALCULATE 




DELTAU- 

— 

TIME LAG BETWEEN POINTS. 




BFREQ - 


FREQUENCY AT WHICH CALCULA- 

TIONS ARE TO BEGIN. 




EFREQ — 

— 

FREQUENCY AT WHICH CALCULA- 

TIONS SHOULD BE STOPPED. 

INPUT 

CARDS 

© 

SKIP-- 



FREQUENCY THAT IIS TO BE 

OMITTED FROM CALCULATIONS. 




M 



IF EQUAL TO 0, IT MEANS 


TO READ ANOTHER SKIP 
VALUE; IF EQUAL TO 1, THEN 
ALL SKIP VALUES HAVE BEEN 
READ IN AND THE SKIP LOOP 
IS TERMINATED. 



TABLE 5-con‘t 


INPUT DATA © 


FREQ 

GAIN 

PHASE 


THIS DATA IS THE 
OUTPUT FROM THE 
IN APPENDIX I. 


OUTPUTS 


PRINTED: 

NLAGS DELTAU BFREQ EFREQ 

XXXX X.XXXX XX.XXXX XX.XXXX 


SKIP FREQ 
X.XXXXXX 

i 

i 

I 

X.XXXXXX 


(LSKIP) 


XX 

i 

i 

i 

XX 


DIRECT 

PROGRAM 



o 

o 


CO 


Q_ 


O 


m 


o 

LjJ 

cn 

Li- 

CD 


UJ 

_J 

CD 

C 

h- 


CO 

CO 

UJ 


>- 

o 

UJ 


o 

UJ 

cc 


UJ 

CO 

<c 

3C 

Q_ 


<t 

CD 


>- 

o 

z 

UJ 


a 

UJ 

cr 




1 

1 




1 

1 


X 

X 









X 

X 



1 




1 


X 

X 









X 

X 








X 

X 

X 



1 




1 

>- 

X 

X 



1 




1 

an 











X 

X 



1 






X 

X 



1 




1 


X 

X 



1 




1 


X 

X 

X 

X 

1 

1 


X 

X 

1 

1 


X 

X 

X 

X 



X 

X 



X 

X 

X 

X 

1 


X 

X 

1 


X 

X 

X 

X 


UJ 

X 

X 


— > 

X 

X 

X 

X 


GO 

X 

X 


<£ 

X 

X 

X 

X 

1 

<x. 

□c 

X 

X 

1 

h- 

X 

X 

X 

X 


Q- 

X 

X 


1 

X 

X 

X 

X 



X 

X 



X 

X 

X 

X 

1 


X 

X 

1 


X 

X 

X 

X 

1 


X 

X 

1 


X 

X 


X 

X 

1 


X 

X 

1 


X 

X 

X 

X 



X 

X 



X 

X 

X 

X 

1 


X 

X 

1 


X 

X 

X 

X 


z 

X 

X 



X 

X 

X 

X 



X 

X 


>- 

X 

X 

X 

X 

1 

<x 

X 

X 

1 

X 

X 

X 




CD 


• 

1 

cc 



X 

X 



X 

X 



X 

X 

X 

X 

1 


X 

X 

t 


X 

X 

X 

X 

1 


X 

X 

1 


X 

X 

X 

X 

1 


X 

X 

1 


X 

X 


X 

X 

1 


X 

X 

1 


X 

X 

X 

X 


>” 

X 

X 



X 

X 

X 

X 

1 

o 

X 

X 

1 


X 

X 

X 

X 


z 

X 

X 



X 

X 

X 

X 


UJ 

X 

X 


ZD 

X 

X 

X 

X 

1 

ZD 

X 

X 

1 

<£ 

X 

X 



1 

o 



1 

I— 



X 

X 


UJ 

X 

X 



X 

X 

X 

X 

1 

or 

X 

X 

1 


X 

X 

X 

X 

1 

u_ 

X 

X 

1 


X 

X 

X 

X 

1 

1 


X 

X 

1 

1 


X 

X 
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TABLE 5 - OUTPUTS con't 
PUNCHED OUTPUTS 

PUNCHED FOR SPECIAL PLOTTING PURPOSES 


TAU 

12 

RXY 

13 

-TAU 

14 

RYX 

xx.xxxx 

1 

1 

1 

1 

1 

1 

1 

INTEGER 

EITHER 

-1 FOR 

NEGATIVE 

NUMBER, 

OR 0 FOR 

POSITIVE 

xx.xxxx 

1 

1 

1 

1 

1 

1 

1 

SAME 

AS 

12 

xx.xxxx 

1 

1 

1 

1 

1 

1 

1 

SAME 

AS 

12 

xx.xxxx 

1 

1 

1 

1 

1 

1 

1 

1 

xx.xxxx 

NUMBER 

xx.xxxx 


xx.xxxx 


xx.xxxx 



TABLE 6 
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FOURIER TRANSPOSES PROGRAM 


DIMENSION C1M0?),C2 14021 , SKIPI50I 


THIS PROGRAM Will TAKE THE FINAL SPECTRA AND SUPPOSSEDLY COMPUTE 
THE CORRELATION CURVE FROM THE SPECTRA 


READ 1, NI.AGS.DELTAU.RFRFO.EFREO 

1 FORMAT I I 4 . 3F1 0 . 0 1 
PRINT 31 

31 FORMAT [4X,5HNLAGS>5X,6HDELTAU,5X,5HBERFQ,5X»5HFFRE0] 
PRINT 37 , NL AGS , DEI TAIJ . BERER, EFREQ 

32 FORMAT [4X,F5.3,5X,F6.4,5X»F5.3,5X,F5.3,///1 

L S K I P = 0 

PRINT 33 ...... 

33 FORMAT f7X,10HSKIP FRFQ . 1 

DO 2 LM.50 , ’ ' 

LSK I P = l. SK I P + 1 

READ 3, SK I P 1 1. 1 > M 

3 FORMAT f F 1 2 . fl . 6 6 X » 121 
PRINT 30. LSKIP.SKIPIU 

30 FORMAT I I7.5X.Ffl.61 
IF £ M T 4.2,4 

2 CONTINUE 

4 CONTINUE 

DO 5 1=1, ML AGS 
Cl fll =0 . 0 
02 t I 1 = 0.0 

5 CONTINUE 
, PRINT 34 

34 FORMAT f ///, IX . 21H FRFQ LESS THAN 8FRFQI 
DO 17 .1=1 , Nl AGS 

READ 6. FREQ, GAIN, PHASE 

6 FORMAT MX,F0.n,1AX,F8.O,8X,FB.n] 

PRINT 76, FREQ, GAIN, PHASE 

76 FORMAT fiX.3F1E.6I 

IF [ FREO-BFRER I 17,7,12 
1? CONTINUE 

7 CONTINUE 
PRINT 35 

35 FORMAT [///,57X,4HFREO,t1.X,4HGAIN.9X,5HPHASEI 
DO 11 K = 1 , Ml AGS 

READ 6, FREQ. GAIN, PHASE 
PRINT 77, FRED, GAIN, PHASE 

77 FORMAT r 5 0 X ,:3F 1 5 . 6 I 

IF f FREO-EFREO ! 8,13,8 

8 CONT1MUF 

DO 9 1=1 , 1 SKIP 
IF [ FRFQ-SK IP II. I ! 9,111,9 
ill. PRINT 78, FREQ, GAIN, PHASE 
28 FORMAT t IX , 25HFREO GA IN_PHASE__I EFT 0UT.3F15.6] 

GO TO 11 

9 CONTINUE 
T AU = n . n 



... TAUN = ‘TAU , - 

DO in 1 = 1 , NJL AGS 

RF.= I2,i .*3 ^L415.9_* Efi E J3 *J.A.U1 , , ..... . . 

RM=f 2. *3 .1 4159*FRFQ*TAUN1 

A=RP-PHASE ..... : . 

B = RM-PH.4SF 

CP = G A I W*COS t A } ... . ........ 

CN=GAIN*COSIB] 

Cl U l=.C.l.C..l l.tGJL 

C^E I 1=C2 ( ? 1+CN 

tau=t au+dfuta u . . 

TAtJN = -T AIJ 

I 0 CONTINUE: . ..... 

II CONTINUE 

13 T A U5.0.J). 

T A U N = * T A U 

.PRINT 41 

41 FORNAT [14X,3HTAU»14X»3HRXY»13X.4H-TAU»14X»3MRYX1 
DO 24 .1=1 , NL APS 

PRINT 40. TAU.C1 m ,TAUN,C?tIl 

40 format !5x.4Fi5.6] , , 

if [ci r i n ' 14, i«5,i5 

14 ci.m=ci m*t-i. j 

i?=-i 
GO TO-. 16 

15 I 2 = fl 

16 IF [TAI.INJ. 17,18,1.8 • 

1 7 T A(IN = T A U N * r -1 . I 

I 3 = - 1 
GO TO 19 

18 13*0 

19 IF 1021111 20,21,21 
?n c? 1 1 1 =C2 m*r-u 

l 4 = - 1 
GO TO 22 

21 14=0 

22 PUNCH 23, T At) , I 2 , 01 f I] » I 3 , T AUN , I 4 , C2 [ IJ 

23 FORMAT f 5X , F.l 0 . 6 , 3 f 4X . I 2, FI 0 . 6 I] 

tau = tau+ofitaij 
taun=-tau 

24 CONTI NUF 
STOP 
END 
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APPENDIX III 


PROGRAM FOR ACCUMULATIVE STANDARD DEVIATION 


This program is designed to calculate the standard deviation of 
piecewise spectral data using the accumulative last piece of data as 
the mean. Other calculations include a standard error calculation for 
each piece of data, and a student's t-distr ibution of the spectral data. 

While the program itself is not complex, it has been included so 
that anyone wishing to use it will not have to write it. 
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1 C. STD - COMPUTED. FROM PIECEWISE-SPECTRA-CaLOH^A-TJONS-BY .USIN 6 -.-EJ.NaL— — 

2 c ACCUMULATIVE SPECTRA 

3 C * *#**##*»**##*# »»»#***#*» * «**»###»* — *.« a 

4 C aaaaaaas##*******#*** 

_-5 DIMENSI0N„-G Ift023 j CQ t»t)2X^QI»O^I-»S Tl> G-t40 g - ) #S-TOe44023^STTeO-tA02J 

6 READ 1 # N/MAX#DELFRQ 

_7- 1 — FORMAT — tl-X i -l 3#2X*I3*3X# £A*OJ — 

8 D 8 20 I * 1 * MAX 

— 9 BEAD-40/-CP-tI4*Q C l ] » G M 3 — — 

10 10 FORMAT [9X#3F8.53 

- 1 1 STDCtn - 0.0 — — 

12 STDQU3-0.0 

.13 STDG Cl J «Q.O 

14 20 CONTINUE 

45 YaMAX-1 i — — 

16 FREQsO.O 

4. 7- —PR I NT— 43 ; ______ 

18 43 FORMAT [5x# 10H SC #4X#10H SG #4x#10H SQ 3 

43 00-50-4*4*^1 

20 SC-0.0 

24 SQ * 0 • 0 — . — — — — — 

22 SG-0‘0 

.23 DO-30 — 1*1* MAX — — — — 

24 READ 40/ COSP# QUAD# GA f N 

25 40— FORMAT— C&X# 3F8 .53 1 

26 STDC C 13 - CCOSP-CO Cl ] 3 **2+STDC (13 

-2-7 STDQCU- CQUAO-Q Cl43*»2+STOCm — 

2 d STDGII3- tGAIN-etIn**2+STD0m 

.23 SC?SC+ CCOSP-CO CI3 3 **2 : 

30 SQ«SQ+tQUAD-Qtm#*2 

.41 SG *S G+-CG A 1 N "G'C133-**2 — 

32 30 CONTINUE 

33 SC sSQRT F CSC/Y3 

34 SQ-SQRTFCSQ/Y] 

-35 SG*SQRTF CSG/Y3 ■— — — 

36 PRINT 41# SC# SQ# SG 

37 PUNCH-44*- SC*SQ*SG 

38 41 FORMAT ClX# 3C4X#Fl0*63 3 

-39-v- 50-CONT I-NUE- 

40 PRINT 110 

-44 — — 44O-F0RMAT- -C5X# 7H FREQ #4X/7H.CQ #4X# 7H -TC #4X#'7H STDC-#-4-X* 

42 17H Q #4X#7H TQ #4X#7H STQQ #4X#7H G #4X#7H TG 

-43 24X#-7H — S-TDG—3 : 

44 DO 60 1 * 1 # MAX 

-45— — X = N — — 

46 STDCU3-SQRTFtSTDCCl]/CX33_ 

,47 STDQ C 1 3 *SQRTFTST-Dq U 3y-CX3-J 

48 STDG Cl] -S qRTF CSTDG CI3 / CX3 3 

-49 TC 3 C O-t-4 3 7STDC-C-H 

50 TQ-QCI3/STDQCI3 

-54 — TG-G C l ) /STDG C 1-3 — 

52 PRINT 111 / FREQ#C 0 cn#TC#STDCtI]#Qtn#T 8 #STDQ[! 3 #GCn#TG#STDG t n 

-53 — -- m -FORMAT-C-lX#- 1-0C4X/-F7 » 4 3-3 
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>0 Co Si Ov 


. „ PUNCH 112^ FREQ* CO [13 *TC*STDC H3 *Q [IiJ *T6 #STdQ £f J *G [II *TG#STDG [13 

U2...F0 RH At...-UX*-FA8-g J-9f 8 * 3> 3 I I __JL 

FREQsFREQ+DEL,FRG 

-6(l-C0HT-lHUE . - 

STOP 

— — END 
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APPENDIX IV 


PROGRAM FOR MEAN STANDARD DEVIATION 


This program is designed to calculate the standard deviation of 
piecewise spectral data by using mean values calculated in the program. 
Also included is a standard error calculation for each piece of data , 
and a student's t-distribution of the mean values to its standard 
deviation. 

This program is included as a time-saving device for those who 
wish to use it. It is essentially the same program that is included 
in Appendix III. 
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1 C STD COMPUTED FROM. P.IECEWISE SPECTRA CALCU^A^flNS_.U^tMCL--M£AKL 

2 C VALUES CALCULATED IN' THIS PROGRAM 

3 C * # * * * **** *# # # a.*.#*.#.*.*.#..*-#. il*_* #«« # *_#.# # »#.** #»*-**j*.s*»- 

4 C *#*# *♦####*#*#######*#■*#* «■*#### *»#**■*#**#.# ###**#•»»***##**##*#### 

— 5 D IMENSICN G [ 4023/ C0X4-024/Q44021/STDG 1 4 021/S-TQC-C 4 023-/ 5T D 044023 

6 Read 1/ N/MAX/DELprg 

— 7 Uf-ORMAt— Cl X » 1 3/3X/4^3/-3X/ P<M>3 — — — 

8 DO 33 I. 1/MAX 

—9 STOCU3«a.O . - 

10 STD3U3-0.0 

-11 ST0QfU«0«O 

12 cecn*o.o 

-13 — CCiJaO.O — 

14 Gcn« 0.0 

-45 33— CONTINUE — — : — — — 

16 Z-N 

-47 - FREQ <0.0 ; 

13 PRINT 43 

-49 43-E0RMAT_c5x^40N SC /-4X/-10 H SC AX/-10H SO- 3 

20 DO 50 J.l/N 

-2 1 SC 1 0 • 0 — 

22 SQ«0»0 

_ 2_3 SG « 0* 0 : 

24 DO 30 I« 1 /MAX 

-25- READ— AO/— COSP/QUAp/GA-tN 

26 40 FORMAT [9X/3F8.53 

-27 — COXI 3" COL I3+C0SP — _ — — _ 

28 OCn»QCn+QUAD 

-29 G C-I1 s G £ I 3 +G A I N — - — : ___ 

30 STDCCI]»STDCCn+C0SP**2 

-3-1 STDQ C U «-S-TDG 1 1 3 +QUAD*»2 , ; * : 

32 STDGLU »STDGtn+GAIN «#2 

-33 SG^5C+-CC-0SP « CCO L4J-/ZT--3-* * 2 -- 

34 SQ*SQ+CQUAD-CQU3/Z33**2 

_35 SGsSG+-C-GA-I-N* CG £147-243 **2 — 

36 30 CONTINUE 

-37 Y«-CMAX- 13- 

38 SC»SQRTFCSC/Y] 

-39— SQaSORTF CSQ/Y1 — . 

40 SG=SQRTFCSG/Y3 

-4-1 PRINT— 4-4/— SC/ SG/SG : 

42 PUNCH 41/ SC#SQ>SG 

-4-3 4-1 -FORMAT — [lX/-3L4X/FlG-«4>43 — - 

44 50 CONTINUE 

-45 PRINT 140 ' ^ - — ^ : 

46 110 FORMAT [5X/7H FREQ / 4 X/ 7H XC /4X/7H TC /4X/7H STDC /4x# 

-47 1741 XQ — /-4X/-7W — TQ /4X/7H- - STQG - / 4 X /-7W XQ- - /A X/ 7W TG 

48 24X/7H STDG 3 

49 -D O 6 0— l u l /NAX- 

50 STDCCl3»SQRTFt[STDCCl3/Z3-CCOtl3/Z3**23 

-54 — STDat 1-3 aSQRTF t tST.pQ C X 3 /-Z3 * LQ U 3 /Z3 « *2] 

52 STDGU]«SQRTFUSTDGCl 3 /Z 3 -[GEn/Z 3#*23 

-53 : X C » C P X n/Z , — — 
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54 

-SB- 


56 

-57- 


3/Z 

3/7- 


XQaQ C 
-XG«G i 

TC*CCOiCI3/Z3/STOCCI3 
JQsCQCI3/Z3/STDQ Cl 3 


58 TG» CGCH/ZJ/STOGCn 

-59- PRINT-lil^-FREQ^XC^-T^SinCX-P^^fQ^S-tMai^a^J-S-TOaC-M- 

60 Hi FORMAT UX*iOC4XiF7*4]3 
-64 R-UNCH— 4447#- 


62 112 FORHAt ClX' F4 « 2* 9F8* 4] 

63 EREQiFREQ+DELFRQ 

64 60 CONTINUE 

-65 STOR- 


66 


END 
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APPENDIX V 


COMPUTER USAGE 


This section of the paper is included to inform others of the 
procedure for operating the small computers at MSFC (the GE-205, 

SDS-930, and CDC-3200). 

The GE-205, located in building 4732, is operated by two Northrop 
employees. Anyone wishing to use this computer must give his program 
to one of the operators, who will then run it for him. Those who wish 
to use the computer themselves must obtain special permission from 
Mr. John Heaman. 

There are three SDS-930 computers available for our use. These 
computers, located in buildings 4200, 4202, and 4610, are self -operated . 
There are special control cards which must be placed at the front and 
back of the user's program (see Table 7). Table 8 gives the operation 
procedures for these computers. A copy of these instructions is usually 
found on the computer console. 

The CDC-3200 computer is located on the ground floor of building 
4200. It is operated by CDC personnel only. NASA has only three hours 
of priority each day: 0800-0900, 1200-1300, and 1500-1600. NASA con- 

tractors have priority for the remaining time. To use this computer, 
you must have a job number and submit a job card with each run. This 
computer, like the SDS-930, has certain control cards which are placed 
before and after the program (see Table 9). 

Data for for the program are stacked the same way for all three 
computers, directly behind the last control card. On the 3200 computer 
a card with a seven - eight punch in the first and second columns should 
be placed behind the data. This is just in case the computer has not 
correctly been told when to STOP. 


209 



TABLE 7 



© 



iTABLE 8 


* * * SDS COMPUTER' OPERATING INSTRUCTIONS * * * 

THIS GUIDE to !-‘ A CHINe OPERATIONS APPLIES to the TYPICAL job which 
CONSISTS OP 

(A) COMPILE A FORTRAN SYMBOLIC DECK 

(13) LOAD and execute the program 
OTHER TYPES OF- COMPUTER RUMS may REQUIRE MODIFIED PROCEDURES. 

I o log in. 

ii. check machine set up, normal set up is as follows 

(A) SYSTEM TAPE ON UNIT with IUNIT SELECT! ON ZERO, 'DENSITY 
SELECT' ON 200, UNIT IN 'AUTOMATIC', 'READY LIGHT' ON. 

(B) SCRATCH TAPE ON UNIT WITH 'UNIT SELECT' ON ONE, 'DENSITY 
SELECT' ON 2C0, UNIT IN 'AUTOMATIC', 'READY LIGHT' ON. 

(C) THE LETTFR 'O SHOULD BE DISPLAYED IN THE COMPUTER CO'nISOLE 
'REGISTER' WINDOW. 

(D) CONSOLE 'PJN-IDLF-STEP' SWITCH SHOULD CE IN 'IDLE' POSITION 

(E) PRINTER SHOULD HAVE 'READY' LIGHT ON. 

III. INITIALIZE MACHINE FOR YOUR RUN 

A. LONG METHOD 

— t. '"check system tape- 'load point' light, if light is not on 

PROCEED AS FOLLOWS < 

(A) PRESS 'START' BUTTON 0" COMPUTER CONSOLE. 

(B) PRESS THE SMALL BUTTONS “N .THE COMPUTER CONSOLE 

NUMBERED 7, 11,12, ?0* 

(Cj PRESS .RuN-lDLE-STFP. S--ITC 'STEP' AND BACK TO . 

- - 'IDLE', WAIT FOR Sy?TE v TAP uf)AD POINT' LIGHT 

2. PRESS 'START' BUTTON ON COMPUTER S ' •',* 

3* PRESS BOTH i M EMORY CLEAR' BUT TO! , SIMULTANEOUSLY. 

PRESS 'START' BUTTON AGAIN. 

5. MOVE 'RUN-IOLE-STEPi SWITCH TO 'RUN'. 

6. PRESS 'FILL-MAG TAPE' SWITCH DOWN (THIS SWITCH IS ON LEFT 

SIDE OF CONSOLE). 

7. PUT YOUR DECK IN CARD READER AND PRESS CARD READER ! READY* 
BUTTON. 

g. WHEN 'INPUT' LIGHT ON RIGHT SIDE OF THE TELETYPE KEYBOARD 

COMES ON TYPE IN 

AC*CR. 

— THEN PRFSS 'RETURN' BUTTON ON TELETYPE keyboard. 

9. THIS MESSAGE TELLS THE SYSTEM CONTROL PROGRAM THAT CONTROL 

MESSAGES ARE ON CARDS IN THE CARD READER. YOUR JOB 

PROCESSING WILL NOW START. 

B. SHORT METHOD 

THIS METHOD WORKS MOST OF THE TIME BUT WILL SOMETIMES FAIL DUE 
TO THE WAY THE PRECEEDING JOB WAS TERMINATED. IF THIS METHOD 
FAILS TRY THE LONG METHOD. 

1. PRESS 'START' BUTTON ON COMPUTER CONSOLE* 

2 • PUT 'RUN-IDLE-STEP' SWITCH TO 'STEP' AND THEN To 'RUN'* 

3. PUT YOUR DECK .IN CARO READER AND P t. S CARD READER 'READY* 

BUTTON. ' 

, - A. IF THE PROCESSING of YOilR -JOB DOES 1 *T. START PUT iRUN-IDLE- 

4 STEP ' SWITCH TO 'IDLE' AND tRY LOv eTHOD. 

IVV at END C YOJR JOB PUT * RUN- IPLE-STEP ' SWITCH TO 'IDLE' POSITION. 

V.’, LOG OFF A VS CLEAN UP YOUR CARDS, PAPER, ETC, 



TABLE 9 
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